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SOUND  PROPAGATION  THROUGH  A  TURBULENT  ATMOSPHERE: 

EXPERIMENTAL  TECHNIQUES  AND  DATA  ANALYSIS 

Final  Report 

Henry  E.  Bass,  Lee  N.  Bolen,  and  John  Noble 
PARGUM  Report  #87-02 

Abstract 

This  is  the  first  of  a  series  of  reports  on  a  three  year  study  of  sound  propagation  through  a 
turbulent  atmosphere.  This  report  documents  the  experimental  configuration  and  describes  data 
analysis.  The  data  analysis  includes  plots  of  the  real  and  imaginary  parts  of  the  acoustic  pressure 
as  a  function  of  time  (scatter  plots),  probability  of  observing  a  particular  amplitude,  and  the  more 
familiar  structure  functions. 

A  preliminary  analysis  of  data  suggest  reasonable  agreement  in  structure  functions  at 
frequencies  of  500  Hz  and  above.  At  lower  frequencies,  phase  and  log  amplitude  structure 
functions  are  both  larger  than  predicted  from  theory.  A  tentative  explanation  for  this  difference  is 
under  development  and  will  be  presented  in  the  third  of  the  three  volume  senes.  The  second 
volume  will  be  devoted  to  refractive  effects. 


Sound  Propagation  Through  a  Turbulent  Atmosphere: 

Experimental  Techniques  and  Data  Analysis 

I.  Introduction 

of  sound  waves  close  to  the  ground  is  a  complex  problem  involving 
many  interesting  mechanisms.  In  addition  to  geometrical  spreading  and  molecular  absorption, 
which  are  reasonably  well  understood,  the  three  main  mechanisms  which  influence  the  acoustic 
field  are  reflection  with  phase  change  due  to  the  finite  impedance  of  the  ground,  refraction  by  wind 
and  temperature  gradients,  and  scattering  by  atmospheric  turbulence.  Outdoor  sound  propagation 
in  a  turbulent  medium  is  not  a  well  understood  process  and  has  only  recently  begun  to  receive 
serious  attention.  This  report^ will  present  experimental  data  taken  under  several  different 
meteorological  conditions  including  various  degrees  of  turbulence  and  different  geometrical 
configurations.  This  report  also  contains  the  results  of  a  preliminary  analysis  of  the  data. 

The  first  section  will  describe  the  experimental  configuration  of  the  experiments  which 
were  undertaken  over  a  period  of  two  years.  A  complete  set  of  geometrical  configurations  will  be 
presented  along  with  the  analysis  procedure  for  processing  the  acoustical  data.  A  unique  technique 
for  examining  the  fluctuations  in  acoustic  amplitude  is  described  in  this  section.  This  analysis  was 
completed  for  each  microphone  in  the  array  over  a  range  of  frequencies  starting  at  62.5  Hz  and 
going  up  to  8000  Hz.  The  next  section  describes  the  meteorological  information  that  is  available 
for  the  experiments.  Temperature  and  wind  velocity  vary  strongly  with  height  within  the  first  few 
meters  above  the  ground.  Because  of  this,  the  wind  speed,  wind  direction,  and  temperature 
measurements  are  given  as  a  function  of  height.  Sound  speed  profiles  have  been  computed  for 
each  of  the.  runs  where  weather  information  is  available.  The  third  section  deals  with  extracting 
sound  pressure  levels  from  the  analyses  described  in  Section  H 

The  fourth  item  discussed  will  be  the  form  of  the  complex  amplitude  of  the  sound  wave  as 
it  propagates  through  a  turbulent  medium.  As  an  acoustic  wave  propagates  a  distance  r  through  the 
atmosphere  from  a  point  source  S  to  a  receiver  R,  atmospheric  turbulence  causes  fluctuations  in 
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the  acoustical  refractive  index,  the  effect  of  which  is  to  produce  fluctuations  in  the  phase  and 
amplitude  of  the  sound  field  at  the  receiver.  In  order  to  observe  these  fluctuations,  a  set  of  plots, 
referred  to  as  scatter  plots,  are  made  which  show  the  fluctuation  of  the  complex  amplitude  over 
time.  These  plots  will  provide  clues  about  the  interaction  between  the  acoustic  wave  and  the 
turbulence.  From  these  scatter  plots,  an  empirical  argument  emerges  which  provides  a  statistical 
distribution  haying  the  same  characteristics  as  the  data.  A  comparison  between  the  statistical 
distribution  and  the  data  is  also  presented. 

Finally,  the  phase  and  log-amplitude  structure  functions  are  calculated  for  the-  data  runs. 
The  calculations  are  compared  with  the  structure  functions  determined  by  Daigle. 

n.  Experimental  Configuration 

Measurements  were  made  over  a  relatively  flat  open  farm  land  over  a  period  between  mid- 
June  1984  and  mid-July  1985.  The  sound  source  was  driven  by  a  tape  which  had,  pre-recorded,  a 
signal  consisting  of  a  mixture  of  seven  pure  tones  centered  at  one  octave  spacing  beginning  with 
62.5  Hz.  A  run  consisted  of  an  eight  minute  record  of  signals  received  simultaneously  at  five 
microphones  mounted  one  meter  above  the  ground  surface  and  one  microphone  mounted  near  the 
source.  The  received  signals  were  recorded  on  a  seven  channel  TEAC  R-81  tape  recorder. 
Typically,  at  least  one  run  was  made  in  the  morning  and  one  run  in  the  afternoon.  This  allowed  for 
propagation  through  different  turbulence  conditions.  During  a  typical  run,  the  first  five  channels 
were  used  to  record  the  signals  from  the  array  microphones,  the  sixth  was  used  to  record  the 
reference  microphone  and  channel  seven  was  used  for  a  voice  log. 

The  measurements  were  made  using  four  different  types  of  geometries.  The  most  common 
geometry  had  the  source  on  the  ground  with  a  transverse  array  of  five  microphones  one  meter 
above  the  ground  about  100  meters  from  the  source  and  one  microphone  five  meters  from  the 
source  also  one  meter  above  the  surface.  Another  common  geometry  had  the  source  on  top  of  a 
100  foot  tower  with  a  transverse  array  of  five  microphones  on  the  ground  with  one  microphone 
about  five  meters  from  the  source  at  the  same  height  as  the  source.  One  experiment,  December  13 
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-  Run  3.1,  was  performed  with  the  source  on  top  of  the  tower  and  a  vertical  array  of  microphones 
ranging  from  0  meters  to  2.57  meters  above  the  ground.  Another  run,  December  13  -  Run  2.1, 
was  conducted  with  the  source  on  top  of  the  tower  with  a  longitudinal  array  of  microphones  on  the 
ground.  A  complete  list  of  geometries  for  all  of  the  runs  is  found  in  Appendix  A. 

The  transverse  array  was  chosen  so  that  the  phase  and  log-amplitude  structure  functions 
could  be  calculated  for  various  transverse  separation  distances.  This  gives  a  measure  of  the 
fluctuations  in  the  sound  field  after  propagating  some  distance  along  almost  parallel  paths. 

One  quantity  of  interest  is  the  fluctuation  in  the  amplitude.  In  order  to  observe  these 
fluctuations,  we  devised  the  analysis  procedure  shown  in  Figure  2.1.  The  analysis  begins  by 
retrieving  the  data  from  tape  and  passing  it  through  a  1/3  octave  filter  set  at  one  of  the  broadcast 
frequencies.  The  output  from  the  filter  was  sent  to  a  multichannel  analyzer  (MCA)  which  stored 
the  amplitude  of  each  acoustic  half  cycle  in  a  channel  number  proportional  to  the  amplitude.  The 
number  of  counts  in  each  channel  is  proportional  to  the  probability  that  an  acoustic  half  cycle  will 
have  a  particular  amplitude.  The  MCA  accumulated  data  over  the  time  of  the  run.  After  the 
accumulation  was  completed,  the  MCA  showed  amplitude  distribution.  After  examining  these  for 
various  frequencies,  we  found  that  the  amplitude  distributions  were  of  two  basic  shapes;  a 
Gaussian  or  Rayleigh  as  shown  in  Figure  2.2  and  Figure  2.3.  From  these  distributions,  the 
relative  sound  pressure  levels  were  calculated  by  taking  from  the  distribution  plots  the  difference 
between  the  mode  of  the  array  microphone  and  the  mode  of  the  reference  microphone.  The 
procedure  for  calculating  the  relative  sound  pressure  levels  is  outlined  in  detail  in  Chapter  IV. 
Later,  the  amplitude  distribution  plots  will  be  used  to  fit  a  statistical  distribution. 

The  next  analyses  yield  the  phase  and  log-amplitude  structure  functions  and  scatter  plots. 
The  configuration  of  the  analysis  equipment  is  given  in  Figure  2.4.  The  output  from  each  channel 
of  the  recorder  (other  than  seven)  was  passed  through  a  1/3  octave  filter  set  at  one  of  the  broadcast 
frequencies  to  a  multichannel  ADF12F  A/D  converter  connected  to  a  Masscomp  5535  mini-main 
frame  computer.  The  data  was  stored  in  a  16  bit  binary  file  format  for  further  analysis.  The 
digitization  rate  was  twelve  times  the  signal  frequency. 
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The  program  phaz.f  in  Appendix  B  is  used  to  calculate  the  phase  and  amplitude  difference 
between  any  two  channels  for  each  acoustic  half  cycle.  It  accomplishes  this  by  reading  in  1024 
points  from  each  channel  at  one  time  and  subtracting  the  dc  bias  from  each  channel.  Next,  the 
program  searches  through  the  data  and  locates  the  sign  changes  on  each  channel.  The  time  of  zero 
crossing  is  then  estimated.  A  similar  procedure  based  on  change  of  slope  gives  the  corresponding 
amplitude  located  between  the  crossings.  After  the  time  for  each  zero  crossing  has  been  calculated, 
the  difference  between  this  time  and  the  corresponding  time  on  another  channel  is  calculated  which 
is  related  to  the  phase  difference  between  the  two  channels.  The  ratio  of  amplitudes  for  the  two 
channels  being  compared  is  also  computed  for  each  half  cycle.  This  process  continues  until  all 
possible  comparisons  between  channels  have  completed.  The  output  from  this  program  is  used  by 
the  scatter  plot  and  structure  function  routines  described  later. 

HI.  Weather  Profiles 

Meteorological  effects  can  have  a  significant  effect  on  the  received  sound  field.  Large 
eddies  are  formed  in  the  atmosphere  by  instabilities  in  the  thermal  and  viscous  boundary  layers  at 
the  surface  of  the  ground.  Futher  instabilities  cause  these  eddies  to  break  down  progressively  into 
smaller  and  smaller  sizes  until  the  energy  is  finally  dissipated  by  viscosity  in  eddies  approximately 
1  mm  in  size.  A  statistical  distribution  of  eddies,  which  we  call  turbulence,  is  therefore  present  in 
the  atmosphere  at  all  times.  The  intensity  of  the  turbulence,  however,  is  strongly  dependent  on 
meteorological  conditions. 

While  the  effects  of  wind  and  temperature  gradients  appear  similar  in  the  graphs  in 
Appendix  C,  the  following  differences  should  be  noted.  Because  temperature  is  a  scalar  quantity, 
the  refraction  of  sound  produced  by  lapse  or  inversion  conditions  is  the  snme  in  all  horizontal 
directions.  Wind,  however,  produces  refraction  which  is  nonuniform  in  direction  according  to  the 
vector  component  of  wind  relative  to  the  direction  of  propagation.  Thus,  the  refraction  produced 


by  wind  is  zero  when  the  sound  propagates  directly  crosswind,  and  increases  progressively  as  the 
direction  of  propagation  deviates  from  this  condition. 

Temperature  and  wind  velocity  vary  strongly  with  height  within  the  first  few  meters  above 
the  ground.  These  vertical  gradients  produce  steep  sound  speed  profiles  close  to  the  ground.  In 
addition  to  the  vertical  gradients,  the  temperature  and  wind  velocity  fluctuate  about  their  mean.  The 
resulting  random  fluctuations  of  the  refractive  index  scatter  sound,  which  leads  to  random 
fluctuations  in  the  phase  and  amplitude  of  the  sound  wave.  Even  when  the  turbulence  is 
sufficiently  weak  that  it  has  negligible  effect  on  the  sound  field  in  free  space  it  is  still  sufficient  to 
affect  the  sound  field  above  the  boundary,  especially  in  regions  of  destructive  interference  where 
the  sound  level  is  critically  dependent  on  phase  relationships. 

For  most  of  the  experiment  measurements  of  the  wind  velocity  and  temperature  were  made 
at  four  heights  simultaneously.  These  points  were  located  at  3, 10, 30,  and  1 10  feet.  An  example 
of  the  output  from  the  sensors  is  given  in  Table  3.1. 


Time 

Height 

feet 

Wind  Speed  (m/s) 
Avg  Max  Min 

Wind  Direction 

Avg  Max  Min 

Temperature  (C) 
Avg  Max  Min 

11:25 

110 

S.O 

11.0 

6.1 

34.8 

404.3 

386.4 

2.7 

3.1 

2.3 

30 

6.2 

7.8 

4.7 

52.3 

423.3 

402.2 

2.9 

3.3 

2.3 

10 

5:8 

7.3 

4.1 

226.5 

248.1 

196.9 

3.1 

3.6 

2.0 

3 

4.4 

6.2 

2.4 

86.0 

457.1 

438.7 

4.2 

4.6 

3.6 

11:30 

110 

7.5 

9.8 

5.5 

35.9 

404.3 

388.0 

2.7 

3.3 

2.1 

30 

6.1 

7.7 

4.9 

54.4 

420.7 

404.3 

2.9 

3.3 

2.5 

• 

10 

5:6 

6.8 

4.0 

213.8 

227.0 

178.9 

3.1 

3.5 

2.7 

3 

4.3 

5.8 

3.2 

87.1 

462.4 

435.5 

4.2 

4.6 

3.7 

11:35 

110 

7.9 

9.4 

6.1 

37.0 

409.1 

384.8 

2.7 

3.3 

2.2 

30 

6.1 

7.8 

4.4 

54.4 

426.0 

404.9 

2.8 

3.5 

2.1 

10 

5.7 

7.8 

3.7 

183.7 

203.2 

176.8 

3.1 

3.5 

2.3 

3 

4.5 

7.3 

2.9 

90.3 

459.2 

439.7 

4.3 

4.9 

3.6 

Table  3.1 
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This  type  data  is  available  for  seven  of  the  experiments.  The  four  points  at  each  height  represent  an 
average. 

The  sound  speed,  sv,  was  calculated  from 

sv  =  331.45  *  Vfl  +  T/273.15J  +  w  *  cos(<)>  -  90)  (1) 

where  T  is  the  temperature  in  Celcius,  w  is  the  wmd  speed,  and  <J>  is  the  angle  that  the  wind  is 
blowing  relative  to  teh  direction  of  propagation.  A  computer  program  performed  a  least  square 
curve  fit  of  the  sound  speed  data  to  an  equation  of  the  form 

sv  =  svO  +  m  *  ln(z)  (2) 

where  svO  is  the  sound  speed  at  one  meter,  m  is  the  slope  of  the  curve,  and  z  is  the  height  above 
the  ground.  Typical  values  generated  by  this  curve  fit  are  given  in  Table  3.2. 


Date 

Run 

Slope 

cO 

December  13, 1984 

1.1 

0.627 

338.347 

1.2 

0.663 

338.609 

2.1 

0.664 

338.185 

4.1 

0.948 

338.770 

January  1 1, 1985 

2.1 

-0.4876 

325.320 

2.2 

-0.4183 

325.394 

Table  3.2 


A  positive  slope  indicates  that  the  atmosphere  is  downward  refracting  and  a  negative  slope  is 
upward  refracting. 

Also,  in  Appendix  C  there  are  weather  data  taken  near  the  ground  with  a  vertical  array  of 
thermocouples  and  temperature  transducers.  For  the  Sandusky  site  and  June  23,  1985  Bondville 
site,  both  types  of  measurements  were  made.  The  measurements  made  with  the  thermocouple  are 
indicated  with  the  word  trun  and  the  other  is  run.  Some  measurements  made  at  the  Flatville  site 
used  only  the  thermocouple  array. 


IV.  Relative  Sound  Pressure  Level 


The  relative  sound  pressure  level  was  calculated  for  each  run.  In  all  the  measurements,  a 
microphone  placed  close  to  the  source  served  as  a  reference. 

The  relative  sound  pressure  levels  were  computed  from  the  distributions  recorded  by  the 
MCA.  The  MCA  output  had  the  form  of  Figure  4.1.  There  are  two  distributions  contained  in  the 
figure.  The  distribution  marked  B  is  the  distribution  of  amplitude  at  an  array  microphone  and  D  is 
the  amplitude  distribution  at  the  reference  microphone.  The  peaks  marked  A,  C,  and  E  are 
calibration  points  arsed  to  determine  the  dB  levels  of  peaks  B  and  D  from  a  calibration  tone  put  on 
each  channel  prior  to  every  run.  The  difference  in  dB  between  the  mode  of  the  amplitude 


distribution  recorded  for  the  array  and  reference  microphone  is  the  relative  sound  pressure  level. 
The  difference  is  calculated  as 


RdB  =  RodB 


2oij!E£^l.2oiog|Ml 

LK6'xrij  [givscij 


where  R0dB  is  the  difference  in  dB  between  the  reference  calibration  point  and  the  array 
microphone's  calculation  point,  xrg  is  the  MCA  channel  number  for  the  reference,  xq  is  the  MCA 
channel  number  of  the  mode  for  the  ith  array  microphone,  xcg  is  the  MCA  channel  number  for  the 
calibration  point  for  the  ith  array  microphone,  gC6  is  the  gain  of  the  amplifier  for  the  reference 
channel,  gq  is  the  gain  for  the  ith  array  microphone  channel,  grg  is  the  gain  for  reference,  and  gq 
is  the  gain  for  the  ith  array  microphone. 

The  following  example  will  show  how  the  equation  is  used.  Refering  to  Figure  4.2,  peak 

A  is  xq,  peak  B  is  xq,  peak  D  is  xcg,  and  peak  E  is  xrg.  Reading  from  the  graph 

xq  =  2.6 
xq  =  6.1 
xcg=  3.2 
and 

xrg  =  5.05. 

From  the  data  logs  ,  gcg  =  grg  and.  gq  =  gq,  also  RodB  =  44-84  =  -40  dB.  Substituting  these 


Lrw-rcxr*  irtjrr x ak -tjuu.-v  "w  \aa\  f.'r^vx-jcrvfVJi  -LW  r. *  >*  k  ■  w  *  v  ^  t  xy  x  u‘>  uh  ltw  ut*  inrin 


a 

y 


values  into  Eq.  (4.1),  is  -42.28  dB.  For  most  runs,  the  calibration  gain  is  the  same  as  for  the 
run  gain;  however,  there  are  a  few  runs  where  the  gains  are  not  the  same.  Plots  of  relative  sound 
pressure  level  for  every  run  are  included  as  Appendix  D. 

V.  Distribution  Functions 

Large  eddies  are  formed  in  the  atmosphere  by  instabilities  in  the  thermal  and  viscous 
boundary  layers  near  the  ground.  Further  instability  causes  these  eddies  to  break  down  into 
progressively  smaller  sizes  until  the  energy  is  finally  dissipated  by  viscosity  in  very  small  eddies. 
A  statistical  distribution  of  eddies  of  various  sizes  is  therefore  present  in  the  atmosphere  at  all  times 
The  fluctuations  in  the  sound  pressure  level  of  a  pure  tone  measured  outdoors,  however, 
are  frequently  much  larger  than  would  predict  on  the  basis  of  sound  propagation  in  an  unbounded 
medium.  The  influence  of  turbulence  on  the  sound  field  can  be  large  when  a  boundary  is  present 
because  the  field  above  the  boundary  is  critically  dependent  upon  the  phase  relationships  between 
the  direct  and  reflected  waves  near  a  boundary,  then  fluctuations  in  phase  and  amplitude  are  not 
independent 
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A  polar  representation  of  the  signal  is  shown  in  Figure  5.1. 

Note  that  the  signal  is  composed  of  many  single  path  vectors.  These  single  path  vectors  are  a 
mafhftmafip.al  representation  of  the  various  fluctuations  in  the  phase  and  amplitude  of  the  acoustical 
field  as  it  propagates  through  the  atmosphere.  The  total  amplitude  is  a  complex  analysis 
representation  of  the  complex  signal  amplitude  P  plus  a  complex  random  noise  contribution.  The 
random  noise  term  takes  into  account  the  random  nature  of  the  fluctuations  in  the  index  of 
refraction  of  the  turbules  and  size  distribution  of  the  turbules  in  the  atmosphere  over  time.  Any 
fluctuations  in  the  total  amplitude  will  depend  upon  the  fluctuations  in  the  single  path  phases  and/or- 
amplitudes,  and  the  noised 

fa  order  to  view  the  change  in  phase  and  amplitude  over  time,  a  series  of  scatter  plots  were 
made.  A  scatter  plot  is  a  representation  of  the  changes  in  the  complex  amplitude  over  a  fixed 
amount  of  time.  Samples  of  these  plots  are  provided  in  Figure  5.2  and  Figure  5.3.  The  samples 
are  from  the  January  11,  1985,  run  2.1  data  set 

These  representative  scatter  plots  appear  as  pictured  in  Figure  5.4. 


Figure  5.4 

Assume  that  the  real  and  imaginary  parts  have  a  Gaussian  distribution  about  averages  x  and  y. 
These  assumptions  give  a  Bivariant  Normal  Probability  Function  of  the  form 


Channel  #  4  Frequency  =125  Hz 


Qjnssajcj  XjdutSduij 


Figure  5.3 
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Converting  to  poiar  coordinates, 


P(p)  =  p/cr 


2 

Iflppoto 


(5.1) 


where  p0  =  the  distance  to  the  center  of  the  distribution,  \xo  +  ^°  / 
c  *  the  standard  deviation  of  the  distribution 
p  =  the  distance  to  a  point  in  the  distribution,  (x2  +  y2) 1/2 
Ifpo/CT»  1,  Eq.  5.1  reduces  to 


P(P)  = 


1 

i2 KG 


e-  (p-pf/2o* 


which  has  the  form  of  a  skewed  Gaussian.  If  Pq/ct  «  1,  Eq.  5.1  reduces  to 


Pip).!.-'''1’1 

2 

a 

which  has  the  form  of  a  Rayleigh  distribution. 

The  mode  of  a  distribution  is  the  point  where  the  maximum  probability  occurs.  This  is 
found  by  taking  the  derivative  of  the  distribution  function  and  equating  it  to  zero.  If  this  is  carried 
out  on  Eq.  (5.1),  the  result  is 


2 

X  -XqX 


*iM 

MM 


- 1=0 


where  x  =  p/cr  and  xQ  -  Pq/c. 


Case  1:  x0  »  1  which  is  supposed  to  be  Gaussian  yields 

Ii  (XqX] 

7-) — f  1  as  x  qX  becomes  large 

ioM 

x2-x,jX-1»0 
x  *  x0  which  means  that 

the  mode  of  the  distribution  is  p0. 

This  result  corresponds  to  that  for  a  Gaussian  distribution. 

Case  2:  Xq  « 1  which  is  supposed  to  be  Rayleigh  yields 
X  (XaxI 

— r  -+0  as  x„x  becomes  small 

IoM 

x2  - 1  =  0 

x  =  1  which  means  that  the  mode  of  the  distribution  is  tr.  This  corresponds  to  the 
behavior  of  a  Rayleigh  distribution. 

Equation  5.1  was  fit  to  the  amplitude  probability  data  from  the  MCA  giving  values  for  a 
and  P0.  Comparisons  between  measured  and  computed  distributions  are  given  at  Appendix  E. 
The  program  used  to  perform  a  least  squares  curfit  to  the  data  is  in  Appendix  B. 

VI.  Structure  Functions 

For  straight-line  propagation  (absence  of  refraction)  a  distance  r  through  atmospheric 
turbulence,  the  log-amplitude  and  phase  structure  of  functions  in  a  plane  perpendicular  to  the 
direction  of  propagation  are  defined  as  * 

Dx(r,p)  =  (jx  (r+p)  -  x  (r 

anc 

Ds(r,p)  =([<(>  (r  +  p)- 4>  (r] 

respectively,  where  p  is  the  separation  between  receivers,  x  is  the  log  amplitude,  and  <|>  is  the 
phase.  The  mean  square  log-amplitude  and  phase  fluctuations,  (x^)  =  ^In  A/A0)2)  and  (s->  = 


<(<j>  -  4>q)2),  where  A0  and  <j>0  are  the  amplitude  and  phase  in  the  absence  of  turbulence, 
respectively. 

If  L  is  a  measure  of  the  scale  of  the  turbulence,  then  in  the  case  (kr)^  »  L,  the  theory 
predicts  that  u  is  the  rms  fluctuations  in  the  acoustic  index  of  refraction 

D1(r,pJ  =  2[(x2)-Bx(p)]  {63) 

Ds(r,p)  =  2[<s2>-Bs(p)]  (64) 


where 


<*%<s2>=(f)<»Va. 


u  is  the  rms  fluctuation  in  the  acoustic  index  of  refraction,  k  is  the  propagation  constant,  and  Bx(p) 
and  Bs(p)  are,  respectively,  the  covariances  of  the  log-amplitude  and  phase  fluctuations .5 

Bx(p)  Bs(p)  q>(p/L) 

(x2)  (s2)  p/L  (6.6) 


where 


0(p/L) 


/•p/L  2 

=  J  e^du. 

* 


In  practice  Eq.  (6.5)  agrees  with  experimental  results  for  (s^).  However,  measurements 
show  that  (x2)  « (s^)  and,  in  addition,  the  log-amplitude  fluctuations  quickly  saturate. 
Substituting  Eq.  (6.6)  and  (6.5)  for  x  into  Eq.  6.3  yields 

Dx(r,p)  =  Vi'<u2)k2rL 

L  J.  (6.7) 

The  extraction  of  the  mean  square  amplitude  and  phase  fluctuations  requires  processing  of 
the  recorded  acoustical  signals.  The  analysis  requires  the  assumption  that  the  phase  and  amplitude 
sequences  represent  samples  of  random  processes  where  ensemble  averages  equal  time  averages 
(Taylor's  Frozen  Turbulence  Hypothesis).  Let  i,  j  refer  to  two  microphones  in  a  plane 


perpendicular  to  the  direction  of  propagation.  Then  the  log-amplitude  for  the  ith  microphone  at  the 
nth  time  sample  is 


xin=ln 


A: 


in 


A. 


10/ 


(6.8) 


where  A^  denotes  the  amplitude  and 


1  N 

A{o  =  2rfAjn 

n=l 


(6.9) 


is  the  average  amplitude  over  N  samples.  The  amplitude  structure  function,  Eq.  (6.1),  is  then 
computed  from 


1  £  2 
^x=  M"^(xin"  xjn)  • 
n-1 


(6.10) 


Similarly  the  phase  structure  function,  Eq.  (6.2)  is  obtained  from 


sfittta-v 

i*l  ' 


2  N  N 
IN 

n«l 


-|2 


(6.11) 


Equation  (6. 1 1)  is  written  in  its  more  general  form  where  the  last  term  accounts  for  the  fact  that  in 
experimental  practice  the  mean  phase  difference  of  the  measured  signals  may  not  be  zero. 

There  are  three  programs  which  are  used  for  this  task.  They  are  called  phaz,  phase,  and 
turbl.  The  purpose  of  the  program  phaz  was  discussed  in  section  2.  The  program  phase  takes  the 
output  from  phaz  and  calculates  the  phase  structure  function  using  Eq.  (6.1 1).  The  program  turbl 
is  used  to  perform  the  calculations  needed  to  arrive  at  the  log-amplitude  structure  function. 

The  output  from  phase  and  turbl  was  plotted  against  transverse  separation  of  the 
microphones.  These  plots  can  be  found  in  Appendix  F.  The  solid  line  is  the  theoretical  cur/e 
based  on  Eq.  (6.7)  using  the  values  from  Table  6. 1  for  (u^)  and  L  measured  according  to  the 
method  outlined  by  Johnson,  Raspet,  and  Bobak.^ 


**  w  VJT*  i VMM  u.xw 


§2 

y 


rnizsi 

Run  2.1  Run  2.2 


mm 4 

Run  1.1  Run  1.2  Run  2.1  Run  4.1 


s 

w 


<u2>  1.6  xlO-6  2.56xl0-6  8.8x10-6  11.7  xlO'6  8.41x10*6  15.4x10-6 

L  15  m  0.82  m  15.0  m  15,0  m  15.0  m  2.17  m 


Table  6.1 


References 

1.  G.A.  Daigle,  J.  E.  Piercy,  and  T.F.W.  Embleton,  "Effects  of  Atmospheric  Turbulence  on 
the  Interface  of  Sound  Waves  Near  a  Hard  Boundary,"  J.  Acoust.  Soc.  Am.  £4,  622-630 
(1978). 

2.  Peter  N.  Mikhalevsky,  "Characteristics  of  CW  Signals  Propagated  Under  the  Ice  in  the 
Arctic,"  J.  Acoust  Soc.  Am.  20,  1717-1722  (1981). 

3.  S.  F.  Clifford  and  R.J.  Lataitis,  "Turbulence  Effects  on  Acoustic  Propagation  Over  a 
Smooth  Surface,"  J.  Acoust  Soc.  Am.  22,  1545-1550  (1983). 

4.  M.A.  Johnson,  R.  Raspet  and  M.T.  Bobak,  "A  Turbulence  Model  for  Sound  Propagation 
from  an  Elevated  Source  Above  Level  Ground,"  J.  Acoust.  Soc.  Am.  SI,  638-646  (1987). 

5.  V.N.  Karavinikov,  "Fluctuations  of  Amplitude  and  Phase  in  a  Spherical  Wave,"  Akast. 
Zh.  3,  175-186  (1957). 


I 

I 

1 

8 

m 


I 


22 


n 


iiaHarRKsnnnDanwwanarajmBsioanaOTTOoniTWTwiTirwramjrT^ji^iw  »-.i.*-jjruw-uj 


B. 


D. 


E. 


F. 


APPENDICES 


Geometries 

This  appendix  contains  the  geometrical  configuration  for  all  of  the  experimental  runs 
performed. 


Programs 

This  appendix  contains  the  computer  programs  used  to  analyze  or  compute  the  results 
obtained. 


Weather 

Contains  the  plots  of  the  different  meteorological  parameters  measured  during  the 
experiment 


Relative  Sound  Pressure  Levels 

Contains  the  plots  of  the  relative  sound  pressure  level  for  each  microphone  and  each 
experiment 


Distribution 

Contains  the  MCA  data  that  has  been  compared  with  the  bivariant  normal  probability 
function. 


Structure  Function 

Has  the  comparisons  between  the  structure  functions  calculated  from  the  data  and  Daigle's 
theoretical  structure  function. 
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APPENDIX  A 


Geometrical  configuration  for  all  of  the  experimental  runs  performed. 
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GEOMETRY 
JUNE  19-20,  1984 


TRANSVERSE  DISTANCES 


1/2 

34.0m 

1/3 

34.1 

1/4 

35.5 

1/5 

41.0 

2/3 

0.1 

2/4 

1.5 

2/5 

7.0 

3/4 

1.4 

3/5 

6.9 

4/5 

5.5 

25 
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Bondville,  III. 
Dec.  13,  1984 
Runs  1.1  &  1.2 


GEOMETRY 


Source  and  ref  mic 
(#6) 


30.48  m 


152.24  m 


Transverse  Distances 


5.0  m 
6.1  m 
6.4  m 
43.0  m 
1.1  m 
1.4  m 
38.0  m 
0.3  m 
36.9  m 
36.6  m 


#5  -  36.6  m< 


#4  -  0  m 


#3  -  0.4  m 


1  #2  -  1.5  m 


#  1  -  6.4  m 


I 

[I 
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Bondville  III. 
Dec  13,  1984 
Run  4.1 


GEOMETRY 


Transverse  Distances 


1/2 

1/3 

1/4 

1/5 

2/3 

2/4 

2/5 

3/4 

3/5 

4/5 


5.0  m 
6.1  m 
6.4  m 
43.0  m 
1.1  m 
1.4  m 
38.0  m 
0.3  m 
36.9  m 
36.6  m 
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BONDVILLE,  ILL 
JAN.  11,  1985 
RUN  2.1 


GEOMETRY 


|  speaker  &  reference  mic  (#6) 


30.48  m  vertical 


91.44  m 


All  microphones  except  ref.  1m  above  ground 


Transverse  Distances 


5.5  m 

6.6 
7.0 
27.0 
1.1 

1.5 

21.5 
0.4 
20.4 
20.0 


#5  -  20.0  m 


#4  -  0.0  m 


#3  -  0.4  m 


#2  -  1.5  m 


#1  -  7.0  m 


Ch.  #  Mic.  # 


voice 


o^VjJVtV*' 


mm 


2/4 

2/5 

3/4 


1.5 

21.5 
0.4 


6 

7 


ref 

voice 


1 
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BONDVILLE,  ILL. 
JULY  23,  1985 
RUN  #1 


GEOMETRY 


Transverse 

Distances 

Ch.  # 

Mic.  # 

1/2 

20.0  m 

1 

7 

1/3 

20.4 

2 

5 

1/4 

21.5 

3 

4 

1/5 

27.0 

4 

3 

2/3 

0.4 

5 

2 

2/4 

1.5 

6 

? 

2/5 

7.0 

7 

voice 

3/4 

1.1 

3/5 

6.6 

4/5 

5.5 

32. 
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BONDVILLE,  ILL. 
JULY  23,  1985 
RUN  #2 


GEOMETRY 


speaker  &  reference  mic  (#6) 


30.48  m  vertical 


745  m 


Transverse  Distances 


m 

u 


BONDVILLE,  ILL 
JULY  23,  1985 
RUN  #3 


GEOMETRY 


m 


Transverse  Distances  Ch.  #  Mic.  # 


1/2 

20.0  m 

1 

7 

1/3 

20.4 

2 

5 

1/4 

21.5 

3 

4 

1/5 

27.0 

4 

3 

2/3 

0.4 

5 

2 

2/4 

1.5 

6 

? 

2/5 

7.0 

7 

voice 

3/4 

1.1 

3/5 

6.6 

4/5 

5.5 

S 
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BONDVILLE,  ILL. 
JULY  25,  1985 
RUN  #1 
RUN  #2 


GEOMETRY 


Source  &  Ch#6  (3  m  away) 

t _ 

100  m 


Horn  2  m  above  ground:  Bass  1 .36  m  above  ground 

Source  2  m  above  ground 

All  microphones  1  m  above  ground 


TRANSVERSE  DISTANCES  AKG  MICROPHONES 


ch#/ch# 

distance 

mic# 

ch# 

1/2 

7.0m 

7 

1 

1/3 

6.6 

blank 

2 

1/4 

5.5 

4 

3 

1/5 

27.0 

AKG  standard 

4 

2/3 

0.4 

2 

5 

2/4 

1.5 

Richard's 

6 

2/5 

20.0 
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BONDVILLE,  ILL 
JULY  25,  1985 
RUN  #3 


3 

y 


S 


GEOMETRY 


Horn  2  m  above  ground:  Bass  1 .36  m  above  ground 
Source  2  m  above  ground 
All  microphones  1  m  above  ground 
***  notice  ch#2  and  ch#3  switched  *** 
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TRANSVERSE  DISTANCES  AKG  MICROPHONES 


h#/ch# 

distance 

mic# 

ch# 

1/2 

6.6m 

7 

1 

1/3 

7.0 

blank 

2 

1/4 

5.5 

4 

3 

1/5 

27.0 

AKG 

4 

2/3 

0.4 

2 

5 

2/4 

1.1 

Richard's 

6 

2/5 

20.4 

3/4 

1.5 

3/5 

20.0 

4/5 

21.5 
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APPENDIX  B 


Computer  programs  used  to  analyze  or  compute  the  results  obtained. 


u  u 


this  program  will  read  the  weather  data  and  convert  to  svp 
then  perform  the  statistics 


character#32  nfile 
real  m<3) ,b<3) ,mavg 

di men si  on  ui<4, 10)  ,d<4, 10)  ,  t<4, 10)  ,sv< 4, 10)  ,ad< 4)  ,x <4) 
di men si  on  xs<4> ,x2s<  4) , ys<4) ,y2s<  4> ,xys<  4>,n<4),h<4),r2<3) 
pr i nt», " Input  data  file" 
read*,nf i 1 e 

open< 1 ,f i 1 e*nf i 1 e , status*" ol d" > 

rewind  1 

KM  Sail  Cl 


i  <•  w  miu  a 

h< 1 >»1 1 0 
h<2>*30 
h<3)«10 
h<4)=»3 
x< 1 >“3.5571 
x<  2)*2. 3084 
X<3>*1 . ! 147 
x<4)=*-0. 09431 
do  100  i»l ,3 
nd*0 

ad< i )*0 .0 
do  200  j*l ,4 

read<l,*)  w<  i  , j ) ,d< i , j ) ,  t<  i  ,  j > 
ad< i )*d< i , j)+ad< i ) 
if  <d< i , j ) ,ne .8 .0)  then 
nd=nd+l 
endi  f 
conti nue 

ad< i >»ad< i >/f 1 oat<nd> 
do  480  j=l ,4 

sv< i  ,  j  )  »  331.5  *  sqrt<l  +  t  < i  ,  j  >  /  273.15) 

!  w(i,j)  *  .:o-i('ad(i)  -  98)  /  0/. 290/7?) 

if  < t<  i  , j ) .eq .0 .0 .and.w< i , j )  .eq  .0 .0)  then 
sv< i , j )*0 .0 
endi  f 

write<6,508)  h<j),sv<i,j> 
conti nue 


cont i nue 
format  <5x ,2f8 .2) 


c - - - c 

c  this  section  starts  the  statistical  portion  of  the  c 

c  program,  first  we  calculate  the  sums  for  each  subset,  c 


do  700  i=t ,3 
xs< i )=0 .0 
x2s<  i  )=*0 . 0 
ys<  i  )  — 0 .0 
y2s< i )=0.8 
xys< i )=0 .0 
n< i )=4 


VW rfV «V/V  : a. v.'  jv\.vvv«r r*  :v.v  - j-  a* 
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750 


@00 


do  800  j**l , 4 

if  <sv< i , j) .eq ^0 .0)  go  to  750 
xs< i )»xs< i ) +x<  j  > 
x2s< i >*x2*< i >+x<j)*x< j) 
ys< i )=ys< i >  +sv< i  ,  j  > 
y2s< i >*y2s< i ) +sv  < i ,  j ) #sv  < i , j ) 
xys< i >=xys< i )  +  su ( i , j  >  *x<  j  > 
if  <sv< i , j ) .eq .0 .0)  then 
n< i )*«n< i )  -  1 
endi  f 
con  t i nue 

m< i )*<f 1 oat<n< i ) >*xys< i  >-xs<  i )*ys<  i  > )/ 

<f  lo*fc<n(  i  >)  *x2s<  i  >-x*<  i  >*xs<  I  >> 
b< i )»<ys< i >-m< i )#xs< i ) )/f 1 oat<n< i ) > 

r2<  i)=<f1oat<n<  i  ) )  *b<  i )  »ys<  i >  Kfl oat<  n<  i ) >  :nn<  i  >  vx  /•=.<  i ) 
-ys< i >*y*< i))/<f1oat<n<i>  )*y2s< i )-ys< i >*ys< i ) ) 

cont inue 

if  <n<l> .eq.n<2> .and.n<l> .eq.n<3>)  go  to  1100 
wr i te<6," ("unequal  data  sets,  check  input")") 
go  to  2008 
1100  continue 

c - c 

c  this  section  calculates  the  combined  statistics  c 


708 


l  300 


1500 

2000 


mavg*0 . 0 
bavg*0 .0 
r2avg*0,0 
nm*8 

do  1300  i*l,3 
mavg*mawg+m< i ) 
bavg*bavg+b< i ) 

r2avg=r2aug+r2< i ) 
nm«nm+ 1 

*  oii !;  i  nue 

mavg»mavg/f 1 oat<nm) 

! .  s.-.j  -! ...'jij, '  f  1  «/.*  t  ( nm> 

r2avg*»r2aug/f  1  oat<nm> 
wri te<6, 1500)  maug ,bavo ,r2aug 


f ormat< lx , 

stop 

end 


5mavg=»"  ,f8.4, "  bavg*"  ,f3.4, "  r2avg=" ,f 3.5) 
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C 

C 

c 

C 

C 

c 

c 

c 

c 

c 


c 

c 


******************************************************** 
*************************************  »<'•'■*  ■  *  ■  :  * 
**********  phaz.f 

**********  Finds  the  amplitude  and  phase 
**********  difference  between  two 

**********  different  channels. 

******************************************************** 
******************************************************** 


********** 


program  phaz 


integer  hiwin,  i,  iblk,  chi,  ch2,  dgfreq,  ipts,  j 

integer  k,  lowin,  nl,  n2,  xcycll,  xcycl2,  nk 

integer  npass,  numch 

integer*2  data(1024,2) 

real  -areal,  area2,  biasl,  bias2,  freq 

real  datal<1024),  data2< 1024) ,  hper,  lastpl,  lastp2 

real  noisel,  noise2,  period,  phasinc,  phastot,  pi,  sumtim 

real  tol ,  tsampl ,  timkpl,  timkp2 

real  t ime< 1024,2) ,  phase<2048),  ampl (1024,2) 

parameter  <p i=3. 141592654,  to1=0.15> 


open <20 ,f i 1 e='/data/adc .  dat' ,  accessed  i rect 

open< 1 ,f i 1 e='/data/hol der" , statu s="old') 

open(2,f  i  1  e='/data/amp 1  .dat' ) 

open(4,f i 1 e=Vusr/dat a/phase .dat' ) 

open(3,f  i  1  e= Vusr/data/f  req  . dat' ) 

rewind  1 

read< 1 ,*) , numch 

read< 1 ,*) , dgfreq 

read< 1 ,*) , i p  ts 

read< 1 ,*) ,chl 

read< 1 ,*) ,ch2 


status=/ ol d/ 


rec 


tsampl  *  1  /  fl oat (dgfreq) 
freq  =  float(dgfreq  /  12) 
hper  =1/(2*  freq  ) 
npass  =  ipts  /  (numch  *  1024) 
xcyc 11=0 
xcycl 2=0 
lastpl  =  0. 

1 astp2  =  0 . 
areal  =  0. 
area2  =  0 . 
phastot  =  0. 
iblk  =  0 
nk  =  0 

***********************************  *•*•***..****.* 
**  Read  In  1024  Points  From  Each  Channel  ** 
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c 

c 

10 


c 

c 

c 

c 


15 


c 

c 

c 


29 

c. 


c 

c 

c 


c 

c 

c 

c 


65 


***** 


**  Skip  First  Hundred  Points  ** 

******************************************* 

low  in  »  chi  +  nk  *  numch  +  100  *  numch 

hiwin  =  chi  +  (nk  +  1023)  *  numch  +  100  *  numch 

J  =  1 

do  2  i =1 ow i n , h i w i n , numch 
re ad (20 , rec=i ) ,data( j ,  1  > 
read(20 , rec=i +ch2-chl ) , data( j ,2) 
j  =  j  +  1 
cont i nue 
nk  =  nk  +  1024 

*****  Determine  DC  Bias  (Just  Average  Over 
*****  Complete  Cycles,  1020  instead  o-f  1024) 

********************************************^*H 

bi asl  =  0.0 
bi as2  =  0.0 
do  15  i  =  l ,1020 

biasl  =  biasl  +  i 1 oat(data( i , 1 > ) 
bias2  =  bias2  +  -f  1  oat(data(  i  ,2) ) 
conti nue 

biasl  »  biasl  /  1020. 
bias2  =  bias2  /  1020. 

******************** 

**  Remove  DC  Bias  ** 

******************** 
do  20  i  =  1,  1024 

datal(i)  *  data(i,l)  -  biasl 
data2( i )  =  data(i,2)  -  bias2 
con  t i nue 


nl=xcycl 1+1 
n2=xcycl 2+1 

******************* ******************************** 
*****  Determine  Local  Zeroes  For  Channels  1  &  2  * 

*************************************************** 
call  cycl e(l , dafcal ,nl ,lastpl ,timkpl , areal ,time ,amp1 
+  tsamp 1 ) 

nl=nl-l 

cal  1  cycl e(2,data2,n2, 1 astp2, t imkp2,area2, t ime ,ampl 
+  tsamp 1 ) 

n2=n2-l 

***************************************** 

*****  Zero  Crossings  are  now  in  ***** 

*****  Time  Array.  Compute  Period  ***** 

***************************************** 
sumt im=0 .0 
do  65  i =1 , n 1 

sumt i m=sumt i m+ 1 i me ( i  ,  1 ) 
con  t i nue 


J 
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SB 
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per i od=2 . *sumt im/nl 
f req=l . 0/per i od 

w********************************************************-* 
*****  Check  Signal  to  Noise  Ratio  for  this  Block  ***** 
w****************************************************-***** 
no i se  l=f loattnzerl )/f 1  oat ( n 1 ) *1 00 . 
no i se2=f 1  oat  t nzer2)/f1 oat<n2)*100 . 
if  tnoi sel . 1 t .5. >  then 

if  <noi se2. 1 t .5. )  then 

***************************** 

*****  Compute  Phase  ***** 
***************************** 

phasi nc=<  t i met  1 , 1 ) -t i me  < 1 ,2) )/per i od*360 . 
phase  <l>=phastot+phasinc 
k=mi n0<nl ,n2> 
do  70  i =2 , k 

phasi nc=<  t imet i , 1 >  —  t i me  < i ,2)>/peri od*360 . 
phase  < i >=phase  < i-l)+phasinc 
cont i nue 
do  75  i  =  l  ,  k 

if  < i bl k . eq .0 .and. i . eq . 1 )  go  to  75 
wr i te(2,*)amp1 <  i  ,  1  > ,  ampl <  i  ,2> 
wr i  te < 4 ,*) phase  < i ) 
con  t i nue 

phastot=phase<k) 
el  se 

print*,'Ch  2.  Too  Noisy.  Skip  Data.  Hold  Phase.1' 

endi  f 
el  se 

print*, 'Chi.  Too  Noisy.  Skip  Data.  Hold  Phase.' 

endi  f 

i bl k=i bl k+1 

if  ( i bl k .ge .npass)  go  to  95 
xcycl 1  =  0 
xcycl 2=0 
if  <n2. gt.nl)  then 

********************************** 

*****  More  Cycles  on  Ch  2  ***** 

********************************** 
xcycl 2=n2-n 1 
do  35  i=l , xcycl 2 

t ime  < i  , 2>  =  t i me  <  n 1  +  i ,2) 
ampl ( i ,2)=ampl <nl+ i ,2> 
conti nue 
end  i  f 

if  <nl .gt .n2>  then 

********************************** 

*****  More  Cycles  on  Ch  1  ***** 

********************************** 
xcyc 1 l=n 1 -n2 
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do  90  i =1 , xcyc 1 1 

t i me< i ,1 >  =  t ime  <  n2+ i , i ) 
amp  1 < i , 1 )=ampl <  n2+ i , 1 ) 
continue 
endi  f 
go  to  10 
wr i te<3,*)  freq 
stop 
end 


************************************* 
**  Finds  the  Zero  Crossing  Between  ** 
**  < 0  , y 1 )  And  <dx,y2)  ** 

************************************* 
real  -function  zerox<yl  ,y2,dx) 
real  yl  ,  y2,  slope,  dx 
slope  =  <y2  -  yl)  /  dx 
zerox  =  -  yl  /  slope 
•  return 
end 


subroutine  c yc 1e<nch, da ta, ncyc, 1  as tpt, timekp, area, 

time,  ampl,  tsampl ) 

***************************************************** 
**  Processes  Cycles  in  Each  Block  o f  1024  Points  ** 
**  Records  Time  and  Amplitude  of  each  1/2  cycle  ** 


a**************#********************************-***-** 


i n teger 
real 
real 
real 


i,  ncyc,  nch 

ampl < 1024,2) ,  t ime< 1024,2) ,  data(l@24) 

aftime,  zerox,  thispt,  lastpt,  prod,  timekp,  tsampl 

p i ,  area 


parameter  <p i=3. 141592654) 


do  10  i=l ,1024 

th i sp  t=f 1  oat ( data< i > ) 
prod  =  th  i  sp  t*l  astp  t 
if  (prod.ge.0)  then 

******************************* 
*****  Same  Side  of  Zero  ***** 


******************************* 


timekp  =  timekp  +  tsampl 

area  =  area  +  tsamp  1  *<.  th  i  sp  t  + 1  astp  t  )/2 

if  ( th i spt .eq .0)  then 

****************************** 

*****  Hit  Zero  Exactly  ***** 
****************************** 
t i me < ncyc , nch )  =  timekp 


'\r*xr*iir»rv/v  wwwimwjRiura 


ampl (ncyc ,nch>  =  < p i *area>/< 2*t i mekp ) 
ncyc  *  ncyc  +  1 
timekp  =  0 
area  =  0 
endi  f 
el  se 

**********************  ************************ 

*****  Crossed  Zero  —  End  of  1/2  Cycle  ***** 

********************************************** 

aftime  *  zerox( 1 astpt , th i spt , tsampl > 

timekp  »  timekp  +  aftime 

area  *  area  +  1 astpt*af t ime/2 

t ime(ncyc ,nch>  =*  timekp 

ampl (ncyc ,nch)  *  < p i *area>/( 2*t i mekp) 

ncyc  *  ncyc  +  1 

********************************************* 
*****  Get  Prepared  For  Next  1/2  Cycle  ***** 
********************************************* 
timekp  *  tsampl  -  aftime 
area  »  th i spt*t imekp/2 
endi  f 

************************************** 

*****  Remember  This  Data  Point  ***** 
************************************** 
lastpt  »  thispt 
conti nue 


return 

end 


***#***•****•*■#**■***•*■#  •a*********************************** 
******************************************************** 
**********  phase,  f  ********** 
**********  Finds  the  average  phase  and  ********** 
**********  the  phase  structure  -function  ********** 
**********  f or  a  pair  o-f  channels.  ********** 
******************************************************** 
******************************************************** 


program  phase 

dimension  tot <1080) , tot2<1000> ,ds<1000> 
i n teger*4  num(1000> 

■format  <  15x,el3.6,12x,el3.6) 

formate  lx ,/,25x , 'Corrected  Results  ' ) 

formatei0x,'Overal 1  Avg  Phase  =' ,el3.6,3x , 'Avg  Structure  - 
f  ormat  ( 25x  ,  '  Raw  Data') 

■format  <  23x  ,  'Very  Raw  Data') 

formate 10x , 'Phase  Structure  Func t i on' , 6x , 'Average  Phase',/ 
3  1 3x , ' <  rad**2) ' , 1 8x , '  <rad)  ') 

format( 10x , 'Total  Phase  Structure  =  ',ei3.6) 
formate 10x , 'Frequency  =  ',fll.6> 
open<4,f i 1 e='/usr/data/phase .dat' ,status='ol d' ) 
opend  ,f  i  1  e='/usr/data/f req  . dat '  , status=' ol  d' ) 
open<2,f i 1 e='/usr/data/ph ' ) 
open <3 , f i 1 e='/usr/data/cph ' ) 
open<5,f i 1 e='/usr/dat a/raw' ) 
rewind  1 
rewind  4 
pi =3. 141592654 
conv=p i/130 . 
read< 1 ,*> ,freq 
ns=l 

wr i te<2,7) 
wr i te<2,3) 
wr i te<5, l ) 
wr i te<5 ,3) 

******************************************************** 
*****  Read  Phase  and  Compute  5  Second  Averages  ***** 
******************************************************** 
j  =0 

kt=20*i f i x(freq) 
do  99  i =1 , 6 

read<4,*) , phase 
>  continue 
i  j=j  +  l 
num< j )=0 
tote j  >=0 .0 
tot2< j>=0 .0 
do  20  i=l ,kt 
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28 


c 

c 

c 

c 

25 

26 


23 


2? 
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read<4,*,end=25) , phase 
t=phase*conv 

if  < i .eq . 1 .and. j .eq . 1 >  then 
t0=t 
el  se 

call  correct i on< t , t6> 

endi  f 

tot< j >=tot< j  >  +  t 
tot2< j >=tot2<j)+t*t 
conti nue 

tot2< j)-tot2< j>/f loat(kt) 
tot<j)=tot<j>/f1oat<kt) 
ds<j)=tot2<j)-tot(j >**2 
wr i te(5,2)ds< j ) , tot< j ) 
num<  j  )=*j 
go  to  10 

******************************************************************** 
*****  Identify  Long  Term  Drifts,  First  Find  Average  Phase,  ***** 

*****  Ignore  Records  Less  Than  5  Sec.  in  Duration  ***** 

******************************************************************** 
j=j-l 

call  average<ds,avg, j > 
call  standard<ds,avg,stan , j ) 
var=stan  *  stan 
pr i nt*,var 

if  <var . 1 t .0 . 1 .or . j .  1  t .6)  go  to  23 
call  search <  tot , tot2,ds,num, j ,stan ,aug> 
call  structLre( tot , tot2,ds, j ) 
go  to  26 
do  29  i=l , j 

wr i te<2,2)ds< i > , tot< i ) 
conti nue 
avg2=0 .0 
avgx=0 .0 
sum=0 . 0 
s i g2=0 . 0 
avgy=0 
do  30  i  =*  1  ,  j 

avgy=avgy+tot< i ) 
sum=sum+tot< i )*f 1 oat<num< i ) ) 
si g2=si g2+f 1 oat(num< • )*num< i ) ) 
avgx=avgx  +  f  1 oat<num< i ) ) 
avg2=avg2+<  tot2< i )-<  tot< i >**2> ) 
con  t i nue 

avgx*avgx/f 1  oat ( j ) 
avg2=avg2/f 1  oat  <  j ) 
avgy=avgy/f 1  oat ( j ) 
sum*sum/f 1 oa  t  <  j ) 
si g2=si g2/f 1 oat<  j >-avgx*avgx 
si ope=<sum-avgx*avgy>/si g2 
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b=avgy-s1  ope*avgx 
wr  i  t  e  <  2 , 4 )  av  gy ,  av  g2 
************************************* 

*****  Now  Correct  -for  Drift  ***** 
************************************* 
do  40  i =1  ,  j 

corr=sl  ope*f 1  oat  <num<  i  ) )  +b 
tot2< i )=tot2< i >-2.*corr*tot< i )+corr*corr 
tot< i )=tot< i >-corr 
conti nue 

************************************** 
*****  Save  Corrected  Results  ***** 
************************************** 
wr i te<3,3> 
wr i te<3,3) 
avgl=0 .0 
avg2=0 .0 
avg3=0 .0 
do  45  i — 1 , j 

ds< i  )*tot2< i ) -tot ( i >**2 
avg3*avg3+ds< i ) 
wr i te(3,2)ds< i ) , tot< i ) 
avgl=avgl+tot < i > 
avg2='avg2+tot2<  i ) 
conti nue 

avgl=avgl/f 1 oat< j  > 

ds< j+1  >=avg2/f  1  oat<  j  )-*avgl*avgl 

avg3=avg3/-f  1  oat<  j  ) 

wri te(3,4)avgl ,avg3 

wr i te(3,9)ds< j+1 ) 

wr i te<3, 16)freq 

stop 

end 


subrout i ne  search <mark ,mark2 ,ds,num, j , stan , avg) 

real *4  mark < 1000) ,mark2< 1000) , ds< 1009) 

integer*4  num<1000) 

pi =3. 141592654 

accept=avg  +  stan 

dif-f=*0. 

i  =  1 

i  4  <  i  .  1  e  .  j )  then 

if  (accept .  1 1  .ds<  i  )  )  then 
do  20  k=i ,  j-1 

mark(k)=mark<k+l ) 
mark2<k)=mark2(k+l ) 
num<k)=:num<k  +  l ) 
ds(k)»ds(k+l ) 
con  t i nue 
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j=j  ~  1 

diff=mark<i)  -  mark(l> 
i  =  i  -  1 
e  1  se 

mark( i >=mark< i )  -  diff 

mark2< i >=ds<  i  )  +  mark<i>  *  mark(i) 

endi  f 
i  =  i  +  1 
go  to  10 
endi  f 
return 
end 


subrout i ne  correct i on < phase , ph i 0) 
pi  =  3.14159265353? 

10  flag*l. 

dphi=*phase  -  pht0 
i  -f  <abs<dph  i  >  .gt  .p  i/2. )  then 
i f  <dph i .eq  .0  . )  then 
signal . 
el  se 

si gn=dph i/abs<dph i ) 

endi  f 

phase=phase  -  sign  *  pi  /  2. 
f  1  ag=0 . 
endi  f 

if  <f1ag.eq.0)  go  to  10 
ph i 0*phase 
return 
end 
c 
c 

subroutine  average<x ,avg,num) 
dimension  x<1000) 
sum=0 . 

do  10  i=l ,num 

sum*=sum  +  x  <  i ) 

10  continue 

avg=  sum  /  float(num) 
return 
end 
c 
c 

subrou t i ne  standard< x , avg , stan , num) 
dimension  x<1000) 
diff=0. 
do  10  i=l,num 

diff=diff  +  <x(i)  -  avg)  *  (x<i)  -  aug) 
10  continue 
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dif-f  =  di  f-f  /  -float(num)  ■ 
stain=sqr  t<  di  f- f  > 
return 
end 
c 
c 

subrout i ne  structure <mark ,mark2,ds, j ) 
real *4  markC 1800) ,mark2< 1880) , ds( 1800) 
do  10  i *  1  ,  j 

ds< i )=mark2< i )  -  mark<i>  *  mark<i) 
18  continue 
return 
end 


49 


.1  w-vwu  *#"w  *nu  V^w*  1  wi  r  wrvi ;  war-'  u*  -ww  ~vj*  -tj»  -\,j 


C  *********************************.**********.*.***:***.***.**,*. 

C  ******************************************,►************* 

c  **********  turbl  ******.**** 

c  **********  Finds  the  ampl  tude  structure  ********** 

c  **********  -for  a  pair  of  channels.  ********** 

c  ******************************************************** 

c  ******************************************************** 

c 
c 

program  turbl 

open <  1  ,  f  i  1  e*/ /data/ amp  1  .  dat '  ,  status®"  ol  d-' ) 
open<2,f i  I  e=*//usr/data/aph/ > 
rewind  1 

101  format < 15x , '2048  Point  Average  =  /  ,el3.8) 

102  format< 10x , 'Ampl i tude  Structure  Function  =  ',el3.6> 
suml»0 .0 

sum2*0 . 0 
n*0 

5  read< 1 , *,end*20> ,  ampl 1 ,  ampl 2 

if  < abs< amp  1 1 ) . 1 t . 5 . 0 . or . abs( amp  11 ) . gt . 5090  . )  go  to  5 
if  <  abstamp  1  2> .  1  t .  5 . 0  .or  .  abs(  amp  1  2)  .  gt .  5000  .  >  go  to  5 
n=r<+l 

suml»suml +  abs< amp  1 1) 
sum2*sum2+abs<ampl 2) 
go  to  5 

20  if  <n.ne.0>  then 

avgl^suml/f 1 oatCn) 
au  g2»sum2/ f 1 oa  t  <  n ) 
el  se 

pr i nt*, "There  are  no  acceptable  amplitude  values" 
go  to  50 
endi  f 
rewind  1 
dx=0 .0 
nn=0 

25  suml=0.0 

n=0 

do  30  i*l ,2043 

read< 1 , * , end=40 ) , amp  1 1 , amp  1 2 

if  < abs< amp  1 1 ) . 1 t . 5 . 0 . or . abs< amp  1 1 ) . gt . 5000  . )  go  to  30 
if  < abs< amp  1  2) .  1  t . 5 . 0 . or . abs( amp  1 2) . gt . 5000  .  ^  go  to  30 
n=n  + 1 

ampl l=abs<ampl 1/avgl ) 
amp  1 2=abs<  amp  1 2/avg2> 
ampl l=al og<amp1 1 ) 
ampl 2=al og<ampl 2) 

suml=suml +  <  amp  1 1 -amp  1 2>  *<  amp  1 1 -amp  1 2) 

30  continue 

if  (n  .ne  .0)  then 

dx=dx  +  suml/f 1  oat ( n  > 
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sum^sumlZ-fl  oat(n> 
wr i te<2, 101 )suml 
nn=nn+ 1 
endi  i 
go  to  25 
i i  <n.ne.0)  then 

dx=<  dx  +  suml/-f  1  oat<n))/-float<nn  +  l) 

e  1  se 

dx=dx/-f  1  oat(nn+l ) 

endi  i 

wr i te<2, 102)dx 

stop 

end 
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***** 


*********************  s************************** 

*****  Theory  .  f  ***** 

*****  Generate.-  data  for  the  theoretical  ***** 
*****  plot  of  the  structure  function.  ***** 

A***************************#************.******** 


program  theory 
real  mu2,L,nu , 1 ambda,k 
character*64  xname, yname 
print*, "Name  of  x  file" 
read* , xname 

print*, "Name  of  y  file" 
read*,yname 
open<l,file  =  xname) 
open<2,fi1e  =  yname) 

pr i nt*, "Propagat i on  Distance  in  Meters" 


read*,r 

pr i nt*, "Mu**2" 
read* ,mu2 

pr i nt*, "Outer  Scale  of  Turblence  in  Meters" 


read*, L 

pr i nt*, "Temperature  at  run  in  Celsius" 
read*,T 

pr i nt*, "Frequency  in  Hz" 
read*,nu 

pr i nt*, " Input  maximum  separation" 
read*,sep 

i h i gh  *  i f i x  <  sep  *  10.) 
i  1  ow  =  1 


pi  «  3 . 141592654 
c8  =  331  .6 

c  =  c0  *  sqrt(l  +  T  / 
1 ambda  —  c  /  nu 
k  =  2  *  p  i  /  l ambda 
k2  =  k  *  k 

do  10  i  =  i 1 ow  ,  i  h  i  gh 
rho  =  f 1 oat< i >  /  < 


call  i ntegrate<rho,L,ph i ) 


sf  =  sqrt<pi)  *<u  )*k  *  r  *  L  *  C  1 


where : 
sf 
2 

<u  > 
k 


Structure  Function 


Fluctuating  Index  of  Refraction 
The  Wave  Number 

Straight-Line  Propagation  Distance 
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c 

c 

c 

c 

c 

c- 


10 


c 

c- 

c 

c 

c- 

c 


10 


L 
ph  i 
rho 


Outer  Scale  of  Turbulence 

Integral  -from  O  to  rho  /  L  o-f  exp<-x*x)  dx 
Spatial  Separation  Perpendicular  to  the 
Direction  o-f  Propagation 


s-f  =  sqr  t  <  p  i  )  * 
if  <sf . 1 e .0)  go 


wr i te<  1  ,  *) 


wr  i  te  <  2 ,  *) 
cont i nue 
stop 
end 


rho 

sf 


mu2 
to 
*  L 


* 

10 


K2  *  r  *  L  *  t  1  -phi 


r  h  o ) 


Performs  a  Numerical  Integration  Using  Trapazoid  Pjj  1 
With  Endpoint  Correction 


subroutine  i n tegrat e < rho , L , ph i ) 
real  L 
i tter=20 
upper=rho 
1 ower=0 

dx=<upper  -  lower)  /  float<itter> 
x=1 ower 

cal  1  equat i on<x ,y0) 
y=7  *  y0 

call  der i vat i ve<x ,yp) 
y=y  +  dx  *  yp 
x=upper 

call  equat i on <x ,y0) 
y=y  +  7  *  y0 
call  derivatively,  yp) 
y=y  -  dx  *  yp 
i  =  1 

x= lower  +  dx 

if  (x.  It.  upper  -  dx  /  2.)  then 
if  c 2* i f i x < i  /  2).eq.i)  then 
call  e  q  u a  t i on  < x , yO ) 
y=y  +  1 4  *  y0 
e  1  se 

call  equat i on < x , y0 ) 
y=y  +  1 6  *  y0 

endi  f 
x=x  +  dx 
i  =  i  +  1 
go  to  10 
endi  f 

phi=  dx  *  y  /  15. 
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return 

end 

c 

c - - 

c  Performs  a  Numerical  Derivative  Using  Difference  c 

c  Tables.  This  is  Used  in  the  Numerical  Integration  c 

c  Routine  For  Endpoint  Correction  c 

c - - 


c 

subroutine  der i vat i ve ( x , yp > 
dimension  del ta(20> ,di ff (20) 
dx  =  x  /  100 

xx=x 

do  110  i =1 ,20 

call  equat i on(xx ,de1 ta< i > > 
xx=xx  +  dx 
110  continue 

do  130  j=l ,13 
do  120  i =1 , 20-j 

del ta( i )— del ta( i+l)  -  delta(i) 

120  continue 

di ff ( j )=del ta( 1 ) 

130  continue 
>'P=0 

do  140  i =1 ,13 

if  (dx.ne.0)  then 

yp=yp  +  (-l.)**(i+l)  *  diff(i)  /  (dx  *  float(i)) 
el  se 
yp=0 

end  i  f 
140  continue 
return 
end 
c 
c 

subroutine  equat i on ( x , y) 

y=exp ( -x  *  x) 

return 

end 
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program  di  str  i  b 

character*64  fnam,  xnam,  ynam 
pr  i  n  t* , "  Root  name  -for  data” 

read* , f nam 

xnam  =  "x"//fnam 

ynam  =  " y"//fnam 

op  e  n  <  1  ,  f  i  1  e=x  n  am ,  s  t  a  t  u  s= "  n  ew "  ) 

open<2,f  i  1  e=ynam , status=" new" ) 

pr i nt*, " Insert  p0" 

read* ,  p0 

pr i nt*, " Insert  sigma" 
read*, si gma 

pr  i  nt*,  "Upper  limit  -for  p" 
read*,pu 

pr  i  nt*,  "Lower  limit  -for  p" 
read*, pi 

pr i n t* , " I ncremen t  for  p" 

read*,p i nc 

x0  =  p0  /  sigma 

xu  =  pu  /  si gma 

x 1  =  p  1  /si gma 

xinc  =  pine  /  sigma 

x  =  x 1  -  xinc 

x  =  x  +  x i nc 

t  =  <x  *  x0>  /  3.75 

if  (t.ge.l)  then 

call  bi ugel <x ,x0 , t ,si gma,P) 
e  1  se 

call  bi ul tl <x ,x0 , t ,si gma,P) 
endi  f 

wr i te< 1 ,*>  x  *  si gma 


wr i te<2,*)  P 

if  <x  .  1 1  .xu)  go  to  10 

stop 

end 


subrou tine  bi vge 1 (x ,x0 , t , si gma , P) 
real  10 

10=0.39394223+0 .01328592/ t+<0 .002255319/ < t*t ))-< 0 . 00 1 575*5/ t; 
+  +<0 ,00916231/t**4)-<0 .02057706/t**5)+<0 .02635537/t**6) 

+  -<0 .01647633/t**7)+<0 ,00392377/t**8> 

10  =  10  /  sqrt<x  *  xQ) 

P=  <<x  *  10)  /  sigma)  *  exp<-((x  -  x0)  *  <x  -  x8))  /  2.; 

return 

end 


subrout i ne  bi vl tl <x ,x@ , t ,si gma,P) 
real  10 
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10=1+3 . 51 5622?*t*t+3 . 0399424* t**4+ 1 . 2067492* t**o 
+  +0 .2659732*t**8+0 .0360768*t*^10  +  O .004531  3*t**12 

P  =  <<x  *  10)  /  sigma)  *  exp<-<x  *  x  +  x0  *  x0)  /  2.) 
return 
end 


a 

m 


(Si 


«a 


Pi 
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**************************************************************** 

* 

*  THIS  PROGRAM  IS  DESIGNED  TO  FIT  THE  DATA  USING 

*  LEAST  SQUARES  TO  A  NORMALIZED  B I VARIANT  DISTRIBUTION. 

* 

******************************************  •s************************* 
program  d  i  st-f  i  t 

parameter  <maxd  =  100,  maxp  =  5,  npts  =  580) 
character*!  yesno 
character*64  device 

real  xdata(maxd) ,  ydata(maxd) ,  par<maxp> 
real  xtheory(npts) ,  ytheory< np ts> 
**************************************************************** 


READ  ALL  THE  DATA  INTO  ARRAYS  AND  KEEP  COUNT 


**************************************************************** 
open  <  10  ,  i  i  1  e=/xdata-f  i  1  e' ) 
open  <  1 1  ,-f  i  1  e='ydata-f  i  1  e" ) 
rewind  10 
rewind  11 

**************************************************************** 


EITHER  ASK  FOR  GUESSES  OR  STORE  THEM 


**************************************************************** 
print  *,"best  guess  -for  p0  and  sigma" 
read  *,  par<l>,  par<2) 
npar  =  2 

***##**.*  .s.**.#.****.*****#.****************************************** 


READ  THE  DATA  POINTS 


*******.*.*****.!*.***.*.***.*****.**.*■*.**.************************* 
i  =  1 

10  read  < 10 ,*,end=20)  xdata< i 'i 
read  < 1 1 , * , end=20 )  ydata< i ) 
i  =  i  +  1 
goto  18 

20  ndata  =  i  -  1 
skip  =  0 . 

s****#**#***.*.***.**.*******#.********* 


NORMALIZE  THE  DATA 


s*#*#*********************************************-**-*******-***** 

call  i  ntegrate(xdata,ydata,ndata.> 


UV\^/Vl^^Aa^WlJWW\AnJ7 owwm 


1  JUTiJUIJC'JCfUlff  M^xnxnxn  xn.rcjcn. X-TK."^- 1« XJl'TDTJW 


LOOP,  GETTING  BETTER  PARAMETERS 


30  call  curve-fit  (  xdat  a ,  ydat  a ,  ndata  ,  par  ,  npar  ,  sum) 

********************************************************H 


PLOT  THE  THEORETICAL  CURVE  AGAINST  THE  DATA 
USING  THE  NEW  PARAMETERS 


a-*********************-************************-****************** 
i f  (sk ip .eq .0  .  >  then 


35  print  *,  "which  device?  screen,  plotter,  or  none" 
read  *, device 

i f  ( dev i ce . ne ." none " )  then 

i  f  < dev  i  ce .ne ."  screen" . and. dev i ce .ne p 1 ot ter " )  then 
pr i nt*,dev i ce , " i s  an  invalid  device" 
go  to  35 
endi  -f 

**************************************************************** 


-  CALCULATE  THE  THEORETICAL  CURVE  - 


************************************************************************** 
call  gtheory  < x theory , /theory , par , np ts) 

call  plotdat  < xdata, ydata , x theory , y theory , dev i ce , ndata , np ts) 
endi  f 
endi  f 

************************************************************-**** 


WRITE  OUT  NEW  VALUES  AND  SAVE  THEM 


***■***■****•*»*•##********************■**'*•******■*****•**•■**•********■*■** 
print  *,  "sum,  p0,  sigma" 

print  *,  sum,  (par(i>,  i=l ,npar> 

*  write  (20,*)  sum,  (par(i>,  i=i,npar> 

***■*************************.*.#***#***•*•***■********•#**■*•*•*********.* 

* 

*  -  EITHER  ASK  NEW  PARAMETERS  OR  USE  OLD  ONES  - 


******•#**#***•#*.*.#***#  *#*##•#■**#•*•■*■***•*•■********#.#■*.*#.**.**.**. 
if  ( sk i p . eq  .  0  .  >  then 

print  *,"Do  you  wish  to  try  again?  (y/n/s)" 

read  *,  yesno 

if  (yesno  .eq.  "s">  then 

print*,  "The  number  of  iterations" 


read*,  skip 
go  to  30 
endi  f 

if  (yesno  . eq 
e  1  se 


" y")  goto  30 
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sk  i  p  = 
go  to 
endi  i 
stop 
end 


skip 

38 


subrout i ne  curve* i  t(xdata,ydata,ndata,par ,npar , sum) 
**************************************************************** 


Fits  theory  to  data  by  a  least  squares  method 
invented  by  Guass. 

xdata,ydata  -  data  points  to  fit(array) 
ndata  -  number  o-f  points  in  array 

par  -  array  o-f  quesses  -for  parameters 

npar  -  number  o-f  -free  parameters 


**************************************************************** 
parameter  <maxp  =  5) 

real  xdata(*>,  ydata<*>,  par<*>,  q<maxp>,  u(maxp) 
real  dyda(maxp),  workl(maxp),  work2(inaxp )  ,  s^maxp^maxp) 
real  theory 

**************************************************************** 


ZERO  MAT  S  AND  ARRAY 


**************************************************************** 
do  18  i  *  l ,maxp 

18  q< i )  *  8 

do  20  i  *  1 ,  maxp*maxp 

20  s<i)  »  0 

**************************************************************** 


STEP' THROUGH  ALL  OF  THE  DATA  POINTS 


**************************************************************** 
weight  3  0 . 
sum  =  0 . 

do  30  i  »  1 , ndata 

xval ue  *  xdata< i ) 

ytheo  *  theory<par ,xva1 ue) 

sum  =  sum  +  (ytheo  -  ydata(i))  **  2 

call  der i v(par ,xval ue ,ytheo,dyda) 

do  30  k  =*  1 ,  npar 

q(k)  =  q(k)  +  dyda(k)  *  (ytheo  -  ydata(i)) 
weight  =  weight  +  dyda(k)  *  ytheo 
do  30  m«l,npar 

*  ******************************** 

*  -  INDEX  ON  S  LIKE  2  DIMENSION 

*  ******************************** 

km  *  npar  *  (m  -  1 )  +  k 
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s(km)  =  s<km)  +  dyda(k)  *  dyda(m) 

39  continue 

weight  =  2 .  *  sum  /  we i gh  t 
do  35  k  =  1 ,  npar 

kk  =  npar  *  <k  -  1)  +  k 
s<kk>  =  s<kk)  +  l.  /  weight 

35  continue 

call  gminv  < s , npar , de t , work  1 ,work2> 
call  gmprd  < s , q , u , npar , npar , 1 ) 

*************************************************************.*** 


* 

*  -  vector  U  NOW  CONTAINS  PARAMETER  ADJUSTMENT 

*  - add  THEM  TO  VECTOR  "PAR" - 

* 


***«***********************•************************.****.■********* 
do  40  i  =  1 , npar 

par  <  i  )  =  par  <  i  )  -  u <  i  > 

40  continue 
sum  =  0 

do  50  i=  1 ,  ndata 

ytheo  =  theory< par , xu al ue > 

sum  =  sum  +  <ydata< i )  -  ytheo)  **  2 

50  continue 
return 
end 

real  -function  theory<  par  ,  xval  ue ) 

************************************  ******************************** 
* 

*  THIS  FUNCTION  CALCULATES  A  YVALUE  FOR  THE 

*  DISTRIBUTION 

* 

*********************************************************  *•**#****■*■**•*■* 
real  par<*) 
sigma  =  par<2) 
x0  =  par < 1 )  /  si gma 
x  =  xualue  /  si gma 
t  *  <x  *  x0)  /  3.75 
if  <  t . ge  .  1  >  then 

•;  all  b  i  u  ge  1  <  x  ,  x  0  ,  t ,  s  i  gma  ,  P  > 
e  1  se 

call  b  i  <j  1  1 1  <  x  ,  x0  ,  t ,  s  i  gma ,  P) 
endi  f 

theory  =  P 
return 
end 
c 
c 

subroutine  b i vge 1 < x , x0 , t , s i gma , P) 
real  I 


■JB 

LI 


S3 


% 


m 
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58 

r-i 


60 


nn.n*A«fUlA«JU>'fw nj»UMUPMuiw.mi  mww  n  u 


1=8 .39894223+0 .01328592/ t+<0 .08225531 9/ < t*t > ) -< 0 . 00 1 57565/ t**3) 
+  +<0 .00916231/t**4)-(0 .02057706/t**5)  +  <0 . 02635537/ t  «  <6; 

+  -<0 .01647633/ t**7)+<0 . 00392377/ t**8) 

I  =  I  /  sqrt<x  *  x0) 

P  =  <<x  *  I)  /  sigma)  *  exp<-<<x  -  x0)  *  (x  -  x0>)  /  2.) 

return 

end 


subroutine  bivltl(x,x0,t,si gma , P) 
real  I 

1=1 +3 . 51 56229*t*t+3. 0899424*t**4+ 1 . 2067492*t**6 

+  +0 .2659732*t**8+0 .0360768*t**l 0+0 .084531 3* t • *1 2 

P  =  <<x  *  I)  /  sigma)  *  exp<-<x  *  x  +  *0  *  x0)  /  2.) 
return 
end 

subroutine  der i m< par ,xval ue , ytheo, dyda) 
**************************•*********************.***.****.**.*.***.**..** 
* 

«  EVALUATE  THE  DERIVATIVES  <DYDA)  OF  THE  THE  THEORY 

*  WITH  RESEPCT  TO  THE  ADJUSTMENT  PARAMETERS. 

* 

**********»»#***#««*********************************************. 
real  par<*),  dyda<*) ,  II,  10,  Iratio 
si gma  =  par<2) 
x  =  xualue  /  sigma 
x0  =  par < 1 )  /  si gma 
Z  =  X  *  X0 

I  rat io  =  1 1 (z)  /  I0(z) 

dydad)  =  ytheo  /  sigma  +  <-x0  +  x  *  Iratio) 
dyda<2)  =  1.  -  <x*x+x0*x0)/2 .  +  z  *  Iratio 
dyda<2)  =  -2.  *  ytheo  /  sigma  *  dyda < 2) 


return 

end 

real  function  I0(z) 
real  10 
t  =  z  /  3.75 


if  (t.ge.l)  then 

pol y=0 .39894223+0 .01 32S592/t+0 .0022553 l?/< t*t) 

+  -0 .00157565/t**3+0 .00916281/t**4-0 .02857706/ t^*5 

+  +0 .02635537/t**6-0 .01647633/t**7+0 .00392377/t**3 

10  =  poly 
e  1  se 

10  =  1 .+3.5156229*t*t+3.0899424*t**4+l .  2067492*t**6 
+  +  0 .2659732*t**3+0 . 0360763*t**l 0+0 . 004531 3*t**l 2 

endi  f 


61 


! 
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return 

end 


real 
real 
t  * 
if  < 


function  Il(z) 
II 

z  /  3 

t .  ge . 
pol  y= 


+ 

+ 


e  1 
+ 


II  = 
se 

pol  y- 


.75 

1)  then 

0. 39394223-0. 03?33024/t-0 .00362013/<  t*t> 

+0 .00163801/t**3-0 .01031555/t**4+0 . 02232967/t**5 
-0 .0239531 2/t**6+0 ,01787654/t**7-0 ,00420059/t**8 
pol  y 


end  i 
retu 
end 


II  = 
f 

rn 


0.5+0 .87390594*t*t+0 . 51 493869*t**4+0 . 1 5034934*t**6 
+0 .02658733*t**3+0 .00301 532*t**l 0+0 .0003241 l*t *+12 
pol  y  *  z 


subroutine  gm  i  nv  (  a ,  n  ,  d ,  1  ,m) 

****************#*********************-s (•****■***+*•***■!*•■*■**•*  *  .i. 
* 

*  SUBROUTINE  MINV 


:  * 


* 

# 

* 

* 

* 

* 

* 

# 

* 

* 


PURPOSE 

INVERT  A  MATRIX 


USAGE 

CALL  MINVCA,N,D,L,M> 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


DESCRIPTION  OF  PARAMETERS 

A  -  INPUT  MATRIX,  DESTROYED  IN  COMPUTATION  AND  REPLACED  BY 
RESULTANT  INVERSE. 

N  -  ORDER  OF  MATRIX  A 
D  -  RESULTANT  DETERMINANT 
L  -  WORK  VECTOR  OF  LENGTH  N 
M  -  WORK  VECTOR  OF  LENGTH  N 


REMARKS 

MATRIX  A  MUST  BE  A  GENERAL  MATRIX 


SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 


* 

* 

* 

* 


METHOD 

THE  STANDARD  GAUSS-.JORDAN  METHOD  IS  USED.  THE  DETERMINANT 
IS  ALSO  CALCULATED.  A  DETERMINANT  OF  ZERO  INDICATES  THAT 
THE  MATRIX  IS  SINGULAR. 
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dimension  '  a  < 1 'i , 1 < 1) , m  < 1 ) 

******************************************* -^S*.********:,;*^***.**.***.*.*.* 

* 

IF  A  DOUBLE  PRECISION  VERSION  OF  THIS  ROUTINE  IS  DESIRED,  THE 
C  IN  COLUMN  1  SHOULD  BE  REMOVED  FROM  THE  DOUBLE  PRECISION 
STATEMENT  WHICH  FOLLOWS. 


double  precision  a , d , b i ga , hoi d , dabs 


* 

* 

* 

* 

* 

* 

# 

* 

* 


THE  C  MUST  ALSO  BE  REMOVED  FROM  DOUBLE  PRECISION  STATEMENTS 
APPEARING  IN  OTHER  ROUTINES  USED  IN  CONJUNCTION  WITH  THIS 
ROUTINE. 


THE  DOUBLE  PRECISION  VERSION  OF  THIS  SUBROUTINE  MUST  ALSO 
CONTAIN  DOUBLE  PRECISION  FORTRAN  FUNCTIONS.  ABS  IN  STATEMENT 
10  MUST  BE  CHANGED  TO  DABS. 


* 

*  SEARCH  FOR  LARGEST  ELEMENT 

* 

d=l  .0 
nk=-n 

do  30  k=l ,n 

nk=nk+n 
1 (k)=k 
m<k)=k 
kk=nk+k 
bi ga=a<kk> 
do  20  j=k ,  n 

i z=n*<  j-1 ) 
do  20  i  =k  ,  n 


10 

i j= i z+ i 

if1;  abs(biga)-  abs<a<ij)>)  15,20 

15 

b i ga=a< i j ) 

20 

1 <k)=i 
m<  k )=j 

continue 

************************#■»**•******•!**■*•*•**•*••***•*■*******•*■*■*■**•*****  +  ■«» 
* 

*  -  INTERCHANGE  ROWS  - 


j  =  1  (.  k ) 


4"  •<»  *►  «•»  +  -4b  ^  -4b  ♦  ^  -#»  4b 


25 


i-f(j-k)  35,35,25 
k i=k-n 
do  30  i=l ,n 

k i =k i +n 
hoi d=-a( k i > 
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*  rvxjutrm  v\x  rvxrvjLrv*  AXAX'ATTunurrurTvi 


j i =k i -k+ j 
a.  <  k  i  )=a<  j  i  ) 

30  a<  j i  )  =ho1 d 

*************#***********************************:***-****** :****** 

* 

*  -  INTERCHANGE  COLUMNS  - 

* 

**************************************************************** 
35  i»m<k) 

i  -f  <  i  — k >  45,45,38 
33  j  p=n*< i -1 ) 

do  40  j=l ,n 

j  k=nk+j 
j i=jp+ j 
hoi d=-a( j  k ) 
a<  jk)**a<  j i ) 

40  a<  j  i  >  =ho1 d 

***********************************************'  ***************** 
* 

*  DIVIDE  COLUMN  BY  MINUS  PI MOT  (MALUF  OP  °IMOT 

*  ELEMENT  IS  CONTAINED  IN  BIGA) 

* 

****************************>************************************ 
45  i-f(biga)  43,48,43 

48  d=0 . 0 

return 

48  do  55  i =1  ,  n 

i  -f  <  i  — k  >  50,55,50 
50  ik=nk+i 

a< i k)=a< i k)/<-bi ga> 

55  continue 

**********  -**************************************************** 
* 

*  -  REDUCE  MATRIX  - 

* 

**************************************************************** 
do  85  i =1 , n 

i k=nk+ i 
hoi d=a< i k  > 
i  j  =  i  —  n 
do  85  j=l , n 

i  j  =  i  j  +  n 

i f ( i -k  >  80 , 85 , 80 

80  i  HF  <  j  -k )  82 , 85 , 82 

82  kj=ij-i+k 

a( i j )=ho1 d*a<  kj  >  +  a< i j ) 

85  continue 

**************************************************************** 

* 

*  -  DIMIDE  ROW  BY  PIMOT  - 


ki 

U 


ft 

si 


i 


ra 

to 


n 


i  > 
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k j=k-n 
do  75  ,j  =  l  ,  n 

kj=kj+n 

i-f(j-k)  70,75,70 

a ( k  j  > =a  < k  j  > /b i ga 

continue 


***************************************************** 


-  PRODUCT  OF  PIVOTS  - 


**************************************************************** 


d=d*bi ga 

*****************************************11 


********* 


-  REPLACE  PIVOT  BY  RECIPROCAL  - 


*****************************************************■:•-**■»■*»:•:*■•:-••■•: 
a<  kk )  =  1 . 0/b i ga 
80  continue 

**************************************************************** 


-  FINAL  ROW  AND  COLUMN  INTERCHANGE  - 


**************************************************************** 

k=n 

-  '100  k=<k-l) 

i-Kk)  150,150,105 
105  i — 1 ( k ) 

if(i-k)  120,120,108 
103  jq=n*<k-n 

j  r=n*< i -1 ) 
do  110  j  =  l ,n 

j  k=  j  q  +  j 
ho1d=a< jk) 

•j  i=Jr  +  j 

a<  j  k  >=-a<  j i ) 

110  a<  j i )  =hol d 

120  j  =m  <  k  > 


i-f(j-k)  100,100,125 


.25  ki=k-n 

do  130  i = 1 , n 

k i=k i +n 
hoi  d=a<  k  i  > 

. j  i  =  k  i  —  k  +  j 
a( k i >=-a<  j i ) 

;  30  a< j i )  =hol d 


go  to  100 
150  return 
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end 


subroutine  gmprd<a,b,r ,n ,m, 1 ) 

a*************************************************************** 


SUBROUTINE  GMPRD 
PURPOSE 

MULTIPLY  TWO  GENERAL  MATRICES  TO  FORM  A  RESULTANT  GENERAL 
MATRIX 

USAGE 

CALL  GMPRD<A,B,R,N,M,L> 

DESCRIPTION  OF  PARAMETERS 

A  -  NAME  OF  FIRST  INPUT  MATRIX 
B  -  NAME  OF  SECOND  INPUT  MATRIX 
R  -  NAME  OF  OUTPUT  MATRIX 
N  -  NUMBER  OF  ROWS  IN  A 

M  -  NUMBER  OF  COLUMNS  IN  A  AND  ROWS  IN  B 
L  -  NUMBER  OF  COLUMNS  IN  B 


ft 
ft 
# 

* 

* 

# 
ft 
* 
ft 
ft 
* 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
ft 
* 
ft 
ft 
ft 
ft 
* 

* 
ft 
ft 
ft 
* 
ft 

*ftftftftftftftftftftftft**ft******ftftftj:.Tftftftftft***ftft*ft*ft*ftft***ftft***ftft**ft*ftft**ft* 

dimension  a< 1 > ,b< 1 )  ,r<  1 ) 
i  r*0 


REMARKS 

ALL  MATRICES  MUST  BE  STORED  AS 
MATRIX  R  CANNOT  BE  IN  THE  SAME 
MATRIX  R  CANNOT  BE  IN  THE  SAME 
NUMBER  OF  COLUMNS  OF  MATRIX  A 
OF  MATRIX  B 


GENERAL  MATRICES 
LOCATION  AS  MATRIX  A 
LOCATION  AS  MATRIX  B 
MUST  BE  EQUAL  TO  NUMBER  OF  ROW 


SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 
NONE 


METHOD 

THE 

AND 


M  BY  L  MATRIX 
THE  RESULT  IS 


B  IS  PREMULTIPLIED  BY  THE  N 
STORED  IN  THE  N  BY  L  MATRIX 


BY 

R. 


M  MATRIX  A 


i  k= 

do 


10  K=1 ,1 

i  k«  i  k  +m 
do  10  j=l  ,n 

i  i  •  -  i  r  + 1 
j i=j-n 
i  b3  i  k 
r< i r >=0 
do  10  i =1  ,m 

j i=j i +n 


□ 

$ 

$s 
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i b= i b+ 1 

10  r(ir)=r(ir)+a(ji)*b(ib) 

return 

end 


subroutine  plotdat(xm,ym,xt,yt,device,ndata,npts) 

********************************************************:i :**+*#*  + 

* 

*  Subroutine  to  Plot  Measured  Data  With  Theory 

* 

*******************************************  ******************■**.* 
include  "/usr/i  ncl  ude/1  i  bmp  .  f " 
integer*4  gl s(si zeof gca) 
character*64  device 
real  xm(*) ,ym(*) ,xt(*) ,yt(*) 
integer  bun(2),  1in(2) 
bun(l>  =  0 
bun(2)  =  66 
1  i  nd  )  =  1 
1 in<2)  =  0 

if  (dev i ce .eq . "screen" )  call  system< "cl eargp" ) 
cal  1  mp i n i t  <  g 1 s) 

call  mp lot srcx( gl s, 1 ,ndata ,0 ,xm , " F" , 1,1, null , nu 1 1 > 
cal  1  mplotsrcy(gls,l , ndata,0 ,ym, "F" ,1  ,1  ,nul 1  ,  nul 1 ) 
call  mplotsrcx(gls,2,npts,0,xt,"F",l,l,null , n u 1 1 > 
cal  1  mplotsrcy(gls,2,npts,0,yt,"F" ,1 ,1 ,  nu  1 1  ,  nul 1 ) 
call  mp  1  otchrs( gl s, , 1 ,nul 1 intray, null  intray) 
call  mpl i nes(gl s,2,bun , 1 i n) 
if  (dev i ce .eq ." screen" )  then 

call  mpdev i ce(gl s, "mcd" ,2,0) 
e  1  se 

call  mpdev i ce(gl s, "hpgl 8" ,2,0) 
endi  f 

call  mplot(gls, 0,0,0) 
cal  1  mpend(gl s) 
re ad*, wait 

if  ( dev i ce . eq screen " )  call  system( " c 1 eargp " ) 

return 

end 


subroutine  gtheory  (xt,yt,par,npts) 


* 

*  -  CALCULATES  THE  THEORETICAL  CURVE  - 

* 


real  xt(*),  yt(*),  lolim,  uplim,  par(*),  incr,  diff 
real  theory 


i rujrujrw.*  * jfcM  vj^vr.vK 


Im**' 


1  00 


print  * 

>  "uppe 

r  1 imi t" 

read  *, 

upl  im 

print  * 

>  11 1  owe 

r  limit" 

read  *, 

1  o  1  i  m 

di-f-f  - 

upl im  - 

1  ol  im 

i  ncr  = 

di-f-f  / 

<  np  ts  -  1 > 

xt<l)  « 

1  ol  i  m 

y  t  <  1 )  = 

theory 

<par  ,xt<  1 ) ) 

do  100 

i  =  2, 

np  ts 

xt<  i  ) 

=  x  t  < i -1 >  + 

yt<  i 

=  theory  <p 

i 

cont  i  n 

ue 

i  ncr 


return 

end 


subrout  i  ne  i  ntegrate(xdata,ydata,ndata) 

X#**#**#*##*#**###*#********#**#***##**##***#****#####****#***** 

* 

*  INTEGRATES  THE  DATA  USING  THE  TRAPOZID  METHOD 

* 

*#***##*##**#*#####*###-X*##*###*#*#*#*#**#####***##*##*##***##X# 

real  xdata<*>,  ydata<*> 
area  =  0 . 

do  10  i  »  1 ,  ndata  -  1 

dx  =  xdata<i+l)  -  xdata< i ) 

aughi  =  <ydata(i+l>  +  ydata<  i  )'>  /  2. 

area  *  area  +  dx  *  aughi 

10  continue 

a*********#**#***#***********#******************************-*-*;-**- 

* 

*  -  NORMALIZE  THE  DATA  NOW  - 

* 

do  20  i  -  1,  ndata 

ydata< i )  =  ydata( i )  /  area 

continue 
return 
end 


20 
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program  stordat 

char ac ter *64  root,xroot,  /root 

pr i n t* , " I npu t  root  -for  data  -file" 

read* ,  r  oot 

xroot  =  "xV/root 

yroot  =  "y"//r oot 

open<  1  ,  f  i  1  e--xr oot ,  status*"  new" ) 

open  <  2 , f i 1 e=yr oot , status*" new" ) 

print*, "Input  number  of  points" 

read*,npoi nts 

print*, "Input  position  of  peak  in  dB" 
read*,pdb 

pr i n t* , " I npu t  position  of  peak  from  scale" 
read*,pscal e 

pr i nt*, " Input  height  of  peak" 
read*, he i ght 

pr i n t* , " I npu t  position  of  reference  peak  in  dB" 

read*,pref db 

seal e  =  1 .  /  height 

Pref  =  10**(prefdb  /  20.) 

do  30  i  =  1 ,  npoints 

pr i nt*, " Input  x",i," ,y",i 
read*,p ,prob 

xdb  =  pdb  +  20.  *  logl0<p  /  pscale) 
x  -  10**(xdb  /  20.)  /  Pref 
y  =  prob  *  seal e 
wr i te< 1 ,*)  x 
wr i t  e  <  2 , * )  y 
con  t i nue 
stop 
end 


20  .) 


v  "  i  II 
X  >  I  j 


logl0<p  /  pscale) 
)  /  Pref 


c - 

c  This  program  creates  the  data  -for  constructing 

c  a  scatter  plot 

c - 

c 

program  scatter 

open < 1 ,  f i 1 e=" phase . dat" > 

open < 2 , f  i  1  e=" amp  1  .dat 11 ) 

open  <  3 , f i  1 e="xf i  1  e " ) 

open(4,f  i  1  e="yf  i  1  e"  ) 

open  <20  ,  f  i  1  e="  limits") 

rewind  i 

rewind  2 

xmax=-l 00000 

xm i n=l 00000 

ymax=-l 00000 

ym i n=l 00000 

conv  =  3.14159/130. 

10  read< 1 ,*,end=20)  phase 

read<2,*,end=20)  ampl,amp2 
i f  < amp  1 . ne . 0 . and . amp 2 . ne . 0 . )  then 


amp  =  abs<ampl  /  amp2) 

c - - 

c  Convert  the  data  from  rectangular  c 

c  to  polar  coordinates  c 

c - - 

x  =  amp  *  cos< phase  *  conv) 
y  =  amp  *  si n< phase  *  conv; 

c - c 

c  Find  the  range  of  the  data  c 

c - - 


xmax  =  amax 1 < x , xmax ) 
ymax  =  amax 1 <  y , ymax ) 
xm i n  =  am i n 1 <  x , xm i n ) 
ym  i  n  =  am  i  n  1  <  y ,  ym  i  n  ) 
wr  i  t  e  <  3 ,  * )  x 
wr  i  te  <  4 , *)  y 
end  i  f 
go  to  10 

20  xmax  =  amaxl <abs<xmax) ,abs<xmi n) ) 
ymax  —  amaxl <abs< ymax) ,abs<ymi n) ) 
xm i n  =  -  xmax 
ym in  =  -  ymax 

wr  i  t e  <  20  ,  >  xmax  ,  xm  i  n  ,  ymax  ,  ym  i  n 
stop 


r».ir»irw  u»« wnci^itir  w  vn v wjx n.« x«xi r^m 
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c  ********************************************************* 

c  ***********************************************  '.<<•■  !  ••  * 

c  **********  nuphazscat .  f  ********** 

c  **********  Finds  the  amplitude  and  phase  ********** 

c  **********  difference  between  two  ********** 

c  **********  different  channel  ■*.  ■..•••-■••.!* 

c  **********  For  use  in  the  scatter  program  ********** 

c  ******************************************************  >::* 

c  ******************************************************** 

c 
c 

program  nuphazscat 

integer  hiwin,  i,  iblk,  chi,  ch2,  dgfreq,  ipts,  j 

integer  k,  lowin,  nl,  n2,  xcycll,  xcyc12,  nk 

integer  npass,  numch 

integer*2  data( 1024,2) 

real  areal,  area2,  biasl,  bias2 

real  datal<1024),  data2<1024),  lastpl,  lastp2 

real  noisel,  noise2,  phasinc,  phastot ,  pi 

real  period,  sumtim,  tsampl ,  timkpl,  timkp2 

real  t ime< 1024,2) ,  phase<2048),  amp  1 (1024,2) 

parameter  (p i=3. 141592654) 

open(  1  ,f  i  1  e='hol  der '  , status5*' ol  d' ) 

open<2,f i 1 e=/ampl .dat' ) 

open<  4 ,f i 1 e-' phase . dat x ) 

open <20 ,f i 1 e^'data' , access=" di rect/ , reel =2) 

rewind  1 

read< 1 ,*) , numch 

read< 1 ,*) , dgfreq 

read< 1 ,*) , ipts 

read< 1 ,*) ,chl 

read( 1 ,*) ,ch2 

tsampl  =  1  /  f 1 oat<dgf req) 
npass  =  ipts  /  (numch  *  1024) 
xeye 11=0 
xcycl2  =  0 
lastpl  =  0. 

1 astp2  =  0 . 
areal  =  0, 
area2  =  0 . 
phastot  =  0 . 
iblk  =  0 
nk  =  0 

c  ******************************************* 

c  **  Read  In  1024  Points  From  Each  Channel  ** 

c  **  Skip  First  Hundred  Points  ** 

c  ******************************************* 

10  lowin  =  chi  +  nk  *  numch  +  100  *  numch 


.iTj\njvim.<J'’.»r.Ajw„v\^w:  .vi^iKflxaioi  xn.v--.JLn  H-Krn- *--  *j\  xn.xn*nTtnTi.nx.«xra*iin»^n*nx.j™n:jijii»xx.>r| 


C 

c 

c 

c 


15 


c 

c 

c 


28 

c 


c 

c 

r 


C 

c 

c 

c 
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hiwin  =  chi  +  <nk  +  1023)  *  numch  +  108  *  numch 

j  =  1 

do  2  i — 1 ow i n , h i w i n , numch 
re ad <20 , rec=i ) , data< j , 1 ) 
read<20 ,rec=i +ch2-chl ) ,data<j ,2) 
j  =  j  +  1 
con  t  i  nue 
nk  =  nk  +  1024 

*******************  ********************************* 
*****  Determine  DC  Bias  <Just  Average  Over  ***** 
*****  Complete  Cycles,  1020  instead  of  1024)  ***** 

**************************************************** 
b i asl  =  8.0 
b i as2  =  0.0 
do  15  i=l ,1020 

biasl  =  biasl  +  i 1  oat < data< i , 1 ) ) 
bias2  =  bias2  +  f 1  oat <data< i  ,2) ) 
con  t i nue 

biasl  =  biasl  /  1820. 
bias2  *  bias2  /  1820. 

***< **************** 

**  Remove  DC  Bias  ** 

******************** 

do  20  i  =  1 ,  1024 

datal(i)  =  data<i,l)  -  biasl 
data2< i )  =  data<i,2)  -  bias2 
con  t i nue 


nl=*xcycl  1  +  1 
n2=xcycl 2+1 

******************************************************* 
*****  Determine  Local  Zeroes  For  Channels  1  &  2  ***** 

******************************************************* 
cal  1  cycl e< 1 ,datal ,nl , 1 astpl , t imkpl , areal , t i me ,ampl  , 

+  tsamp 1 ) 

n l=n 1 -1 

cal  1  eye  1 e  <  2 , data2 , n2 , 1 astp2 , t i mkp2 , area2 , t i me , amp  1  , 

+  tsamp 1 ) 

n2=n2-l 

¥:  &  &  #  #  ih  ^  ^  &  ih  -fcr  *«r  •*:  +:  *r  •£  -fc  -4  -4  -4  4  -*•  4  *•*  4  *4  4  *♦■  ♦ 

*****  Zero  Crossings  are  now  in  ***** 

*****  Time  Array.  Compute  Period  ***** 
**********.****************************  *  ■■  *■ 


sumt i m  =  8.8 
do  65  i  =  1 ,  n 1 

sumt i m  =  sumt i m  +  t i me  < i  , 1 ) 


conti nue 

period  =2.  *  sumtim  /  float<nl) 


***** 


t.**^.  ************ 

Check  Signal  to  Noise  Ratio  for  this  Block  ***** 
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end  i  f 
go  to  10 
stop 
end 


************************************* 

**  Finds  the  Zero  Crossing  Between  ** 

**  <0,yl>  And  <dx,y2>  ** 

************************************* 

real  -function  zerox<yl  ,y2,dx) 

real  yl ,  y2,  slope,  dx 

slope  =  <y2  -  yl)  /  dx 

zerox  =  -  yl  /  slope 

return 

end 


subroutine  eye  1 e  <  nch , da t a , ncyc , 1 astp  t , t i mekp , area , 

+  time,  ampl ,  tsampl ) 

***************************************************** 
**  Processes  Cycles  in  Each  Block  o-f  1024  Points  ** 
**  Records  Time  and  Amplitude  o-f  each  1/2  cycle  ** 

***************************************************** 


integer 

real 

real 

real 

parameter 


i  ,  ncyc ,  nch 

ampl < 1024,2) ,  t  i me <  1 024 , 2)  ,  data<Ki24j 

a-ftime,  zerox,  thispt,  lastpt,  prod,  timekp,  tsampl 

p i ,  area 

<pi=3. 141592654) 


do  10  i =1 ,1024 

th  i  spt=f  1  oat<data<  i  ) ) 
prod  =  th i spt*l astpt 
i -f  <prod.ge.0)  then 

******************************* 

*****  Same  Side  o-f  Zero  ***** 

******************************* 

timekp  =  timekp  +  tsampl 

area  -  area  +  tsampl *< th i spt+1 astpt )/2 

if  < th i spt .eq .0)  then 

****************************** 

*****  Hit  Zero  Exactly  ***** 
****************************** 
t i me < ncyc , nch )  =  timekp 
amp  1 < ncyc , nch )  =  < p i *area) /< 2*t i mekp ) 
ncyc  =  ncyc  +  1 
timekp  =  0 
area  =  0 
endi  f 


'7/JrW 


c 

c 

c 


"inc,»w™iJ-«»»\™™r»ininiir*innni»Mr»UTiinnniii«inLuwir>ir«»Tiinmi 


-Ult  ITWlTjr  UT, 


el  se 

^S****************.}*..**..*.***.***#,******!.**.***.,.*.*:*.* 

*****  Crossed  Zero  —  End  o-f  1/2  Cycle  ***** 

*****************•***•*****■****.***  JUS****  *#**.*# 

aftime  =  zerox( 1 astpt , th i spt , tsampl ) 
timekp  =  timekp  +  aftime 
area  =  area  +  1 astp t*af t i me/2 
t i me < ncyc , nch >  =  timekp 
ampl (ncyc , nch)  =  <p i *area)/(2*t imekp) 
ncyc  =  ncyc  +  1 

C  ****************************************  *  ?  i.  * 

c  *****  Get  Prepared  For  Next  1/2  Cycle  ***** 

c  ********************************************* 

timekp  =  tsampl  -  aftime 
area  =  th i spt*t imekp/2 
endi  f 

c  ************************************** 

c  *****  Remember  This  Data  Point  ***** 

c  ************************************** 

lastpt  =  thispt 
10  continue 
c 

return 

end 
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APPENDIX  C 


Plots  of  the  meteorological  parameters  measured  during  the  experiment. 
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APPENDIX  D 


Plots  of  the  relative  sound  pressure  level  for  each  microphone 

and  each  experiment. 
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APPENDIX  E 


Comparison  of  the  MCA  data  with  the  bivariant  normal  probability  function. 
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APPENDIX  F 


Comparisons  between  the  structure  functions  calculated  from  the  data 
and  Daigle's  theoretical  structure  function. 
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